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SECTION  I 


INTRODUCTION 

In  the  static  analysis  of  continua  by  the  finite  element  method, 
the  stiffness  and  load  matrices  are  obtained  by  the  application  of  the 
principle  of  stationary  potential  energy  within  the  displacement  method. 

In  this  approach,  the  strain  energy  of  an  element  is  expressed  in  terms  of 
some  assumed  displacement  functions.  For  the  monotonic  convergence  to  the 
actual  solution,  it  is  necessary  that  these  functions  should  satisfy  the 
compatibility  conditions  at  the  boundaries  of  the  elements.  The  con¬ 
struction  of  displacement  functions  for  the  interior  of  an  element  which 
will  also  ensure  compatibility  at  the  boundaries,  in  most  cases,  is  rather 
complicated  even  for  elements  with  simple  geometry. 

In  1964,  Pian  [1]  proposed  an  approach  for  the  static  analysis  of 
continua,  by  means  of  which  compatible  stiffness  matrices  could  be  obtained 
without  undue  difficulties.  In  Pian's  approach  [2,3,4],  one  assumes  stres¬ 
ses  within  the  element  and  displacements  on  its  boundary,  hence  the  word 
hybrid.  This  approach  is  extended  for  the  nonlinear  analysis  of  multi¬ 
layered  shell  structures  with  shear  deformation  [5].  It  contains  the  fol¬ 
lowing  new  developments  in  finite  element  formulations: 

(a)  the  foundation  of  the  hybrid  stress  finite  element  method  for 
initial  stress  problems  [6], 

(b)  the  adoption  of  the  above  initial  stress  theory  for  shell-type 
structures,  and 

(c)  the  extension  of  the  resultina  shell  theory  for  the  analysis 
of  laminated  shells  with  arbitrary  stress  distribution  across 
the  thickness  coordinate. 
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Since  the  ultimate  goal  of  this  analysis  is  the  solution  of  the  un¬ 
bonded  contact  of  the  tire  with  a  rigid  frictionless  surface,  flat  trian¬ 
gular  reference  surface  elements  are  used.  The  coupling  of  these  flat 
plate  elements  to  form  a  discretized  shell  reference  surface  is  accomplished 
by  the  nodal  displacements  referred  to  a  fixed  global  reference  system, 
which  is  chosen  to  be  the  base  vectors  of  the  toroidal  reference  surface 
of  pneumatic  tires.  Three  rectilinear  displacements  and  two  rotations  are 
assigned  to  each  node  of  the  triangular  surface  element.  Thus,  each  ele¬ 
ment  is  characterized  by  15  generalized  nodal  displacements.  These  surface 
elements  are  built  up  of  layers  of  different  elastic  properties  and  cord 
angle  orientations  with  respect  to  the  local  reference  frame  of  the  element 
under  consideration. 

Within  the  hybrid  stress  method,  an  equilibrating  stress  field  is 
assumed  over  each  layer  in  terms  of  sinple  polynomials  whose  coefficients 
are  then  related  to  the  aforementioned  generalized  displacements  by  the 
hybrid  stress  variational  principle.  Thus,  a  shell  theory  results. 

An  important  aspect  of  the  finite  element  concept  is  that  the  finite 
elements  may  first  be  considered  to  be  disjoint  for  the  purpose  of  ap¬ 
proximating  a  function  locally  over  an  element.  That  is,  one  can  consider 
an  individual  element  to  be  completely  isolated  from  the  collection  of  all 
elements  and  can  proceed  to  approximate  a  function  over  the  element  in  terms 
of  its  values  at  the  nodes  of  the  element,  independent  of  the  ultimate  loca¬ 
tion  of  the  element  in  the  connected  model  and  independent  of  the  behavior 
of  the  function  in  other  finite  elements.  Thus,  it  is  possible  to  develop 
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a  catalog  of  the  various  finite  elements  in  which  nodal  values  of  the  local 
approximation  are  left  arbitrary. 

Clearly,  in  order  to  stay  within  the  philosophy  of  the  finite  ele¬ 
ment  analysis,  the  aforementioned  equilibrating  layer  stress  fields  must  be 
disjoint  from  layer  to  layer.  Moreover,  this  local  (layer)  stress  approx¬ 
imation  must  be  the  same  over  each  layer,  usually  simple  polynomials  whose 
coefficients  wil 1  be  referred  to  as  stress  coordinates.  Recall  now  that  in 
the  usual  finite  element  method,  the  connectivity  of  disjoint  elements  are 
established  by  relating  the  local  coordinates  to  global  coordinates,  often 
called  as  the  assembling  procedure.  In  this  case,  the  connectivity  of  the 
layers  within  a  finite  element  is  established  by  satisfying  the  interlayer 
equilibrium,  which  is  facilitated  by  the  matching  of  appropriate  stress  coor¬ 
dinates  of  the  two  layers  under  consideration. 

Since,  in  the  hybrid  stress  method,  the  interelement  equilibrium 
need  not  be  satisfied  explicitly,  the  set  of  stress  coordinates  are  inde¬ 
pendent  for  each  element.  Connectivity  is  satisfied  in  the  variational 
sense.  Of  course,  the  price  one  pays  for  the  element-wise  disjoint  stress 
field  is  that,  in  the  construction  of  the  element  stiffness  matrix,  one  must 
eliminate  the  stress  coordinates  in  terms  of  the  generalized  displacements 
in  the  hybrid  stress  functional,  in  order  to  construct  the  stiffness  matrix 
associated  with  the  generalized  nodal  displacements. 

The  extension  to  nonlinear  problems  utilizes  an  intrinsic  local" 
coordinate  frame  which  follows  the  element  during  deformation.  The  asso¬ 
ciated  initial  stress  resultants  and  the  effect  of  the  displacements  on  the 
equilibrium  equations  manifest  themselves  in  the  particular  solution  of  the 
element  equilibrium  equations  in  terms  of  the  stress  resultants. 
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The  initial  or  unloaded  geometry  of  the  carcass  is  represented  by  an 
arbitrary  plane  curve  which,  when  rotated  about  the  wheel  axis,  defines  the 
tire  middle  surface.  The  carcass  thickness  varies  in  the  meridian  direc¬ 
tion  but  is  constant  in  the  circumferential  direction.  The  initial  cord 
angle  varies  with  the  meridian  position  according  to  the  classical  lift 
equation  of  bias  tire  construction  [7]. 

The  tire  carcass  material  is  treated  elastically  as  a  laminated 
anisotropic  composite.  Each  ply  is  considered  as  an  orthotropic  layer  con¬ 
sisting  of  elastic  textile  cords  embedded  in  an  elastic  rubber  matrix.  The 
principle  orthotropic  moduli  for  a  given  ply  are  then  calculated  in  terms 
of  the  elastic  properties  of  the  cord  and  rubber  based  on  the  principle 
of  compatible  deformation, 

The  external  loads  applied  to  the  tire  are  inflation  pressure  and 
inertial  forces  due  to  rotation. 

In  the  numerical  implementation  of  the  aforementioned  disjoint 
stress  field  concept,  certain  software  problems  associated  with  the  construc¬ 
tion  of  the  element  stiffness  matrix  were  experienced.  Namely,  for  many 
layered  shells,  extensive  out-of-core  matrix  manipulations  were  required, 
including  the  inversion  cf  the  stress  field  matrix  for  the  elimination  of 
the  stress  coordinates  in  terms  of  the  generalized  nodal  displacments.  These 
out-of-core  operations  resulted  in  prohibitive  computer  time  requirements, 
even  for  a  linear  problem.  Therefore,  alternate  approaches  were  investi¬ 
gated,  which  included  layer  lumping  schemes  and  the  use  of  displacement- 
type  formulations  with  the  capability  of  incorporating  an  arbitrary  shear 
stress  variation  across  the  shell  thickness. 
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In  Section  II,  a  systematic  development  of  initial  stress  formulations 
is  provided,  which  is  then  specialized  for  the  analysis  of  general  toroidal 
shells  in  Section  III.  In  Sections  IV  and  V,  the  appropriate  numerical 
formulations  are  described,  which  include  the  hybrid  stress  finite  element 
formulation  and  a  displacement-type  finite  element  approach  which  is  based 
on  a  piecewise  parabolic  shear  stress  variation  across  the  tire  thickness. 

The  Fortran  extended  version  of  the  Fortran-Compass  computer  code  is 
written  for  the  CDC-6600  machine  under  the  Scope  3.3  and  Scope  3.4.3  oper¬ 
ating  systems. 


SECTION  II 


BASIC  FORMULATIONS  FOR  INITIAL  STRESS  PROBLEMS 

For  clarity  in  presentation,  the  ideas  will  first  be  developed 
here  within  the  framework  of  flat  plate  theory. 

Each  material  particle  in  the  original  configuration  C-j  is  iden¬ 
tified,  in  general,  by  three  curvilinear  coordinates  (i  =  1,  2,  3). 

The  numerical  values  of  which  define  a  particle  in  C^  define  the  same 
particle  in  every  subsequent  configuration  (also  referred  to  as  convected 
coordinates).  To  describe  the  motion  of  the  body  relative  to  C-j ,  a  fixed 
rectangular  cartesian  coordinate  system  x^  (i  =1,  2,  3)  in  three-dimen¬ 
sional  space  is  also  established. 

In  general,  C^  is  defined  to  be  the  configuration  of  the  body  before 
the  addition  of  the  n-th  increment  of  load;  whereas,  CN+1  is  the  config¬ 
uration  of  the  body  after  the  addition  of  the  n-th  load  increment.  In 
configuration  Cfj,  the  states  of  stress,  strain,  and  deformation  are  pre¬ 
sumed  to  be  known.  During  the  process  of  the  n-th  load  increment,  con¬ 
figuration  C^  is  treated  to  be  in  a  state  of  "initial  stress".  Incremen¬ 
tal  displacements  due  to  the  addition  of  the  n-th  load  increment  are  mea¬ 
sured  from  C^.  In  the  following,  as  a  generic  case,  the  movement  of  the 
body  from  the  reference  state  to  the  deformed  state  C^+1  through  small 
but  finite  increments  in  stresses,  displacements,  and  external  loads  is 
treated. 

The  position  vector  of  a  particle  in  is  denoted  by  r  and  that 
of  the  same  particle  in  is  denoted  by  R.  If  x..  are  the  cartesian 

coordinates  of  the  point  in  and  e^  are  cartesian  bases,  it  follows 
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(2.1) 


The  covariant  base  vectors  tangent  to  £.  lines  in  configuration  CN  are  then 
given  by 


9k 


9x. 

Tie' 


(2.2) 


and  the  covariant  and  coitravariant  metric  tensors  and  contravariant  base 
vectors  in  are  given  by 

9rs  "  9r  9S 


(2.3) 


where  6^  is  unity  when  r  =  t,  and  zero  otherwise. 

For  a  flat  plate,  ^  and  are  taken  as  the  reference  surface 
coordinates  and  as  the  thickness  coordinate. 

In  such  a  coordinate  system,  the  vector  field  of  incremental  dis¬ 
placements  from  CN  to  CN+1  is  measured  in  the  basis  system  of  as 

(2-4) 


which  gives  the  usual  representation  of  Green's  strain  tensor  in  CN  as 
ea8  =  1  <ua,6  +  UB.a  +  “v.c  uv,6  +  %  W,B> 
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where 

e  „  =  membrane  strain, 

aB 

u  „  =  curvature  strain,  and 

aB 

y  =  shear  strain,  and 

a 

to  =  total  rotation. 

a 

In  the  reference  state  CN>  which  is  presumed  to  be  known,  let 
the  initial  stress  resultants  be  defined  by 

"aB  =  {  TaB  d53 
maB  °  ^  53  TctB  d?3 

"S  ■  /„  Ta3  d53  (? 


where 

n°„  =  membrane  stress  resultant. 

aB 

m°„  =  couple  stress  resultant, 

aB 

r°  =  shear  stress  resultant, 

a 

=  Piola-Ki rchoff  stresses,  measured  per  unit  area  in  CN,  and 
aB  a3  IN 

h  =  plate  thickness. 

Let  the  initial  surface  loads  and  edge  tractions  in  the  current 
reference  state  CN  be  given  by 
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f°  =  (f°,  p°) 


n  =  (na,  r) 


m  =  (m  ) 

or 


(2.7) 


One  can  then  prescribe  additional  surface  forces  (f  ,  p ) ,  additional 
edge  tractions  (ft  ,  r,  m  )  on  a  portion  S,  of  the  boundary,  and  additional 

LX  LX  I 

displacements  and  rotations  (u  ,  w,  w  )  on  a  portion  S9  of  the  boundary. 

Let  the  corresponding  increment  in  the  Piola-Kirchoff  resultants  be  rep¬ 
resented  by  nag,  mag,  and  ra.  The  principle  of  virtual  work  then  states, 

l  [Ke  +  na6)6eaB  +  <raaB  +  "W6ka8  +  K  +  ra>6\  '  <fa  *  fa)Sua 


-  (p»  +  p)6w]dA  -  /  [<n°  +  iia)«ua  ♦  (m°  +  ma)6a,a 
51 


+  (r°a  +  ra)6w]dS  =  0 


(2.8) 


where  the  strains  are  given  by  Equations  (2.5).  In  Equation  (2.8),  the 
area  A  and  the  boundary  refer  to  the  known  current  reference  state 

and  6  denotes  variations.  Note  that  u  ,  w,  and  u>  vanish  on  S0. 

a  a  2 

Equation  (2.8)  may  be  written  as 

n° 

/  [n  0  6e  Q  +  m  0  5k  0  +  r  5y  +  -Sp-  5(w  w  0  +  u  u  „) 

^  aB  a8  aB  aB  au'a  2  v  ,a  ,B  v,a  v,B' 


+  5(u  a)  0  +  u  0  u)  )  +  r°  6(ui  u  )  -  f  6u 

2  y,ot  y,B  Y,B  Y»a  a  '  y  y  a'  a  a 


pSw]dA  -  /  [na  6ua  +  ma  «u ia  *  r«w]dS  =  -  /  l-f-  S(ua>(J 

^  M 


9 


ju.ii  m.m*  nip  mugii  U  ■  ^wmmmmmrnw 1  ■ 1 .  -m  j.u 


r«wr*i— 


III  Q 

+  UB,a}  +  -f  5(a3a,B  +  V}  +  "a  6(wa  +  W,a} 


-  f°  6u  -  p°  6w]dA  +  /  [n°  6u  +  m°  6co 

a  a  ^  J  i  L  a  a  a  a 


+  r°  6w]dS 


(2.9) 


It  is  now  assumed  that  the  known  initial  stress  state  (n°^,  m°g, 

r°,  p°,  n°,  m°,  r°)  in  CKI  is  in  equilibrium  prior  to  the  addition  of 

a*  a  v  a’  a*  '  N 

the  incremental  loads.  Then,  the  right  hand  side  of  Equation  (2.9)  can 

be  shown  to  be  identically  equal  to  zero.  That  is. 


"Je.e  *  fa  -  0 


m°  q  a  -  r°  =  0 
aB,$  a 


r°  +  p°  =  0 

a, a 


and 


n°  =  n°0  nQ 
a  aB  B 


m°  =  m°D  n0 
a  aB  B 


r°  =  r°  n 

a  a 


(2.10) 


(2.11) 


where  n  (a  *  1,  2)  are  the  components  of  the  edge  normal.  However,  due 


to  the  numerical  incremental  solution  technique  for  solving  the  above 
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<*********** 


fi  ■  inirirY  i ( ■ 1  T|-ti^ffYrr r 


problem,  the  initial  stress  state  may  not  be  in  equilibrium.  Following 
[a],  it  is  possible  to  derive  an  equilibrium  error  check  if  the  right-hand 
side  terms  in  Equation  (2.9)  are  retained. 

Assuming  now  that  the  reference  state  is  one  of  equilibrium, 
Equation  (2.9)  can  be  shown  to  yield  the  equilibrium  equations  for  the 
incremental  resultants  referred  to  the  current  known  reference  state  C„ 
as  follows, 


Ve +  [(\e  *  n;e)ua.6],r  +  [(\e +  "V"«.e].v 

+  [(ry  +  +  fa  =  0 

ma6.B  '  r«  +  C(mY8  +  ^“cAb  *  [(rY  +  ^Vy3 


r  +  p  +  [(n  +  n°Jw  ]  =  0 

a.a  a  8  a8  »a  ,8 


(2.12) 


and 


on  S-| . 

The  principle  of  virtual  work  as  given  by  Equation  (2.9)  can  now  be 
generalized  through  the  usual  methods  into  a  counterpart  of  the  Hu-Washizu 
principle  in  linear  elasticity  [9]. 
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Assuming  that  the  elastic  stress-strain  relations  are  of  the  type 


"aB  =  "au'V  kYX’  YX> 


mctB  =  maB(eY*’  kyX’  V 


ra  =  ra(eYX'  kYX'  YX> 


(2.14) 


or 


ectB  =  eoB(nYX’  "yX*  rX> 


kct6  =  kaB(V  V  rX> 


Y, 


a 


■  W  V  rX> 


(2.15) 


one  can  define  an  elastic  strain  energy  function 


5E  ■  "cB  iec,B  +  “aB  5kaS  +  ra  6yc 


(2.16) 


In  addition,  introduce  the  constraint  conditions  (2.15),  and 


u  =  u 

a  a 


w  =  w 


(0=0) 
a  a 


(2.17) 


on  Sr,. 


Then,  one  can  formulate  the  generalized  functional 

"g  =  "g^aB*  W  V  V  W’  V  naB’  "aB*  ra;  V  ma*  r)  (2,18) 
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where 


\  *  l  (E<ea6'  kc.B>  Yo)  +  -f-  <%  \B  +  Uv,a  Vs> 


in  which 


+  TT  (uY,a  “y,B  +  “y,B  “y.o1  +  ra('Y  Vet1  '  fa  ua 

-  pw  ■  na8^ectB  '  1  (ua,B  +  UB,c.  +  uv,t<  uv,B  +  w,a  ".B1 

« 

•  maB[kaB  '  1  (“a,B  +  V  +  VaV  +  uy,B  “y,«>] 


r  Ty  -  (u  +  w  +u)U  ) ] }dA  -  {  (n  u  +  m  to 
aLra  v  a  ,a  y  y,a'J  i  v  a  a  a  a 

bl 


+  rw)dS  - 


/  [n  (u 

i  L  a  a 
b2 


G  )  +  m  (oa 

a  a  a 


uia)  +  r(w  -  w)]dS 


-  e 


(2.19) 


e  =  -  f  [n°  u  „  +  m°0  to  0  +  r°  (u>  +  w  )  -  f°  u  -  p°w]dA 

b  i  L  aB  a, 6  a8  a, 8  a  '  a  ,<r  a  a  r 


+  J  [n°  u  +  m®  to  +  r°w]dS  (2.20) 

i  L  a  a  a  a 

bl 

In  Equation  (2.10),  e*  is  the  :orrection  tern  to  check  the 
equilibrium  of  initial  stress  state  in  the  reference  state  C^.  If  the 
reference  state  is  theoretically  in  equilibrium,  then  the  variation  of 

e*  with  respect  to  u  ,  to  ,  and  w  is  zero.  It  follows  that  the  Euler  equa- 

ot  ot 

tions  corresponding  to  6Hg  =  0  are: 
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(a)  the  equilibrium  equations, 

(b)  strain-displacement  relations, 

(c)  displacement  boundary  conditions, 

(d)  traction  boundary  conditions,  and 

(e)  stress-strain  relations: 


n 


3E 


aB  3e 


m 


aB 

3E 


'aB  "  IE 


aB 


,  =  9L 

a  9y 


a 


(2.21) 


The  assumption  of  Equation  (2.16)  is  equivalent  to  the  existence  of 
a  potential  B  such  that  % 


B<naS-  "W  ra>  -  naB  eaS  *  "cB  kaB  +  ra  \  '  E 


(2.22) 


Now,  using  Equation  (2.22)  in  Equation  (2.19),  one  obtains  the  functional 

■rrD(n  0,  m  0,  r  ;  u  ,  w,  to  )  such  that 

R'  aB  aB  a  a  a7 

n° 

{  {-  B<"aB’  "aB'  ra)  +  TT  (w,a  W,B  +  Uv.a  “v.B* 


m 


+  -Sr-  (u  a)  +uQu)  )  +  r°(a)U  )  -  f  u  -  pw 
2  v  Y»o  Y,B  Y»“;  a  y  Y»a  a  a 


+  -f  K,i  *  uB,a  +  uv,a  uv.S  +  w,a  ".S*  +  “T  (“a,B 


m 


aB 


+  V  *  \,a  “y.B  *  \,B  “Y.c>  +  rA  +  w,a  +  “y  V«,)dA 
‘  I  <"a  ua  +  "a  “a  +  ™>dS  *  l  A  ‘  Sa>  *  "A  '  “a> 


+  r(w  -  w)]dS  -  e 


(2.23) 
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Now,  the  variation  of  with  respect  to  the  stress  resultants  (n^, 
mag,  ra)  and  the  displacements  (ua>  w;  a^)  yields  the  following: 

(a)  the  equilibrium  equations, 

(b)  displacement  boundary  conditions, 

(c)  traction  boundary  conditions,  and 

(d)  displacement-gradient-resultant  relations: 


1 

2 


<Ua,8  + 


Ue,a  + 


u  ,  u  ,  „  + 

v,a  v,8 


w  w 

,ol  ,r 


1 

1 


(w  Q  + 
a,  6 


0)o  + 

8, a 


UY,a 


UY,8  “y.o^ 


+  w 


,a 


+  OJ 

Y 


(2.24) 


If  in  Equation  (2.19)  one  assumes  a  priori  that  the  strain-displace¬ 
ment  relations  (2.5)  and  the  displacement  boundary  conditions  (2.17)  are 
satisfied,  then  one  arrives  at  the  functional 


7T  —  f  [E(u  ,  W,  W  )  + 
p  i  L  v  or  or 


n  0 

^(w 


m 


W  a 

»a  ,8 


uv,ot  uv^  + 


a8 


^UY,ot  “y.6 


uv,8  UY.a>  +  ra<“Y  \.a>  '  fa  ua  '  Pw]dA  '  l  <"a  ua 


+  m  a)  +  rw)dS  -  e 

a  a  ' 


(2.25) 


which  is  the  potential  energy  functional  for  the  initially  stressed  plate. 

At  this  stage,  it  is  desirable  to  construct  a  theory  which  allows 
large  deflections  but  requires  the  angles  of  rotations  to  be  small  com¬ 
pared  to  unity.  For  this  degree  of  accuracy,  it  is  necessary  to  retain 
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in  the  strain-displacement  relations  only  those  nonlinear  terms  which  in¬ 
volve  the  product  of  the  normal  displacement  gradients,  and  neglect  all 

others.  Thus, 


*  1  (ua,B  *  US,a  +  W,a  ",B> 
koB  =  1  (Vs  +  “B.o5 


V  =  U)  +  W  ^ 
'a  a  ,a 


(2.26) 


Then,  the  generalized  functional  tt  of  Equation  (2.19)  takes  the 


*g  ■  kcB’  Ya>  *  ¥  <%  W,B>  -  fa  Ua  -  pw  ■  "^aB 

'  \  (ua,B  +  UB,u  +  W,«  W.B)]  '  mc.B[kaB  '  7  (“a,B  +  “B,-11 

-  rjL\  '  K  +  W,a)]ldA  ‘  l  ("o  “«  +  "«  “a  +  ^ 


■  /  t"a<ua  -  2a>  +  "A  ’  “a>  +  r<W  ’  ")]dS  "  ' 
S2 


(2.27) 


The  Euler  equations  corresponding  to  the  above  functional  are  the 
equilibrium  equations 


naB,B  +  fa  =  0 


maB,B  ’  ru  =  ° 


ka,a  +  As  +  "asKA  +  P  =  ° 


(2.28) 
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I J. JL  ^  . ,  ■■  ■“*  -*m  <J}»  '•'"  v'.-J 


with  the  boundary  conditions 


n  =  n  0  n„ 
a  a$  8 

ma  =  mag  n6 

?  *  rc  "a  +  (na6  +  naB)w,a  n6 


on  S^,  and 


A 


w 


/\ 


w 

O) 


a 


(2.29) 


(2.30) 


on  S2« 

Furthermore,  the  strain-displacement  relations  as  defined  by  (2.26) 
and  the  stress  resultant-strain  relations  are: 


-  _il_ 

'a8  '  Sea8 


m  -  _2I_ 

a 6  “  8k  c 
aB 


r 

a 


(2.31) 


Next,  certain  simplifications  can  be  made  in  order  to  facilitate 
the  above  developments  for  numerical  calculations.  One  can  assume  that 
each  increment  is  such  that  the  incremental  displacements  u  ,  w,  and 
rotations  aja  are  of  order  0(e);  whereas  the  initial  stresses  are  of  order 
0(1). 
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Then,  the  incremental  strain-displacement  relations  are 


ea$  =  7  (ua,3  +  u3,a)  + 


1 


(“«  a  +  J 


'a6  “  7  %,3  T  “3,o 


Y  =  co  +  w 
'a  a  ,a 


Similarly,  for  elastic  materials, 


nc6  “  0(e)l  na6  “  °<>> 


(2.32) 


"W  “  °<E>;  maB  *  °<’> 

ra  *  0(e);  r°  -  0(1)  (2.33) 


Using  Equations  (2.32)  and  (2.33),  the  third  equilibrium  equation 
in  (2.28)  can  be  simplified  as 


ra,a+  (nSe  V),e  +  P  * 0 


(2.34) 


and  the  third  traction  condition  of  Equation  (2.29)  takes  the  form 


r  =  r  n  +  n°0  w  „  n„ 

a  a  a3  ,a  3 


(2.35) 


With  the  above  simplifications,  the  generalized  functional  tt^,  the 
Reissner  functional  tt^,  and  the  potential  energy  functional  TTp  take  the 
forms 
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"g  *  l  tE<eae>  V  V  +  ¥  (w,o.  W,B>  '  fa  ua  '  ‘  "c,8[ec 


-  \  (ua,8  +  u6,a>]  ’  ma8[ka8  ’  1  <“a,B  *  “8,a)]  '  rA 

-  (ua  +  w>0)]}dA  -  /  (n0  ua  +  Sa  *a  *  fw)dS 

b1 

'  /  ^na^Ua  “  UJ  +  maHi  '  +  r^W  '  w^dS  "  e  (2>36) 

S2 


B(not8*  ma6’  ra]  +  "f  {w,a  W,B}  ‘  fa  Ua  '  pW 

+  ¥  <ua,8  +  uS,a>  +  ¥  <“a,B  +  “m>  +  r«(“c  +  W,a>ldA 
'  I  <"a  ua  +  "a  “o.  +  ™)dS  '  [  ["ce(ua  ‘  Ga>  +  mc.(“a  ’  “a* 


+  r(w  -  w)]dS  -  e 


(2.37) 


"p  *  /  tE<V  “•  “a>  +  ¥  V  W,B  '  fa  ua  ‘  ^  ‘  {  <"a  ua 


+  ma  wa  +  r\,)dS  -  e 


(2.38) 


At  this  point,  note  that  the  functional  of  Equation  (2.37)  leads  to 
a  mixed  numerical  formulation  while  the  functional  of  (2.38)  leads  to  the 
displacement  formulation  of  finite  element  analysis. 

If  in  Equation  (2.36)  one  assumes  that  the  linearized  incremental 
equilibrium  equations,  the  boundary  traction  conditions,  and  the  stress 
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resultant-strain  equations  are  satisfied  a  priori ,  then  the  following 
functional  can  be  constructed  by  utilizing  Equation  (2.22), 

"c  ■  -  l  BKb-  ma$*  ra)dA  +  l  (na  “o,  +  “a  “a  +  ™>dS  '  e‘  (2'39) 

which  is  now  suitable  for  the  hybrid  stress  finite  element  formulation 
with  the  following  modifications. 

★ 

First,  note  that  the  functional  e  is  a  constant  with  respect  to 
variations  in  the  resultants  nag,  ma^,  and  ra»  and  hence  does  not  contri¬ 
bute  any  terms  in  the  variational  equation  6rrc  =  0.  However,  an  equili- 

★ 

brium  check  on  the  initial  stresses  can  be  performed  by  retaining  c  in 
Equation  (2.39),  as  in  the  usual  displacement  formulation.  If  the  condi¬ 
tion 


r  =  r  (2.40) 


on  S-|  is  satisfied  a  priori ,  then  one  can  introduce  it  as  a  subsidiary  con¬ 
dition  and  consider 

\  -  -  l  B(V  ma6’  rc)dA  +  |  <"a  +  ma  “a  +  ™)dS 


+  /  ‘  na)ua  +  '  ma)\  +  ^  “  r)wlds  “  e  (2.41) 

S] 


* 

The  variational  equation  <$7rc  =  0  then  leads  to  the  linearized  strain 
displacement  boundary  conditions  on  Sg  (2.30). 
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It  is  emphasized  that  one  assumes  a  priori  that  the  linearized 
equilibrium  equations  are  satisfied  in  terms  of  the  Piola-Kirchoff  re¬ 
sultants,  taken  per  unit  area  in  as 


"aM  +  fa  =  0 


m«M  '  rc  *  0 


ra,a  +  <"aS  W,a>,B  +  P  =  0 


(2.42) 


Following  Pi  an  and  Tong [4],  one  can  modify  the  functional  in  Equa¬ 
tion  (2.41)  for  the  hybrid  stress  formulation  of  the  finite  element  model 
as 


*  M  , 
TT  =  l  { 

c  n=l 


/  B(n 
^n 


a0'  rna6*  'a 


,  rJdA  +  /  (nx  ua  +  m„  +  rw)dS } 


3A_ 


a  a 


* 


e 


(2.43) 


where  M  is  the  number  of  elements,  Ap  is  the  area  of  the  n-th  element, 
and  3An  is  the  interelement  boundary  of  the  n-th  element.  In  the  above 
equation,  ua,  oja,  and  w  are  the  interelement  boundary  displacements  and 
rotations,  and  the  meaning  of  Lagrangian  multipliers  can  be  used  to  satis¬ 
fy  the  interelement  traction  continuity  requirement  on  the  average.  These 
interelement  boundary  displacements  and  rotations  are  prescribed  such  that 
they  inherently  satisfy  the  interelement  compatibility  conditions. 

The  numerical  facilitation  of  Equations  (2.38)  and  (2.43)  for  the 
stress  analysis  of  pneumatic  tires  will  be  discussed  in  Section  IV. 
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SECTION  III 


INITIAL  STRESS  FORMULATIONS  FOR  GENERAL  TOROIDAL  SHELLS 
A.  Geometrical  Considerations 

The  linearized  field  equations  in  the  presence  of  initial  stresses 
are  developed  here  for  toroidal  shells  with  arbitrary  reference  meridian 
profiles,  following  the  basic  ideas  of  Section  II. 

Let  x.  be  rectangular  cartesian  coordinates  with  the  origin  fixed  in 
space,  as  indicated  in  Figure  1.  The  position  vector  of  a  point  P  is 
defined  by 

r  =  xi  e.  (3.1) 

where  are  unit  base  vectors.  Introducing  the  curvilinear  coordinate 
system, 

X1  =  <ro  +  ?  cos  cos  ^1 

x2  =  (ro  +  5  cos  %)  sin  ^l 
x3  =  k  sin 

C  =  ?U2)  (3.2) 

the  base  vectors  of  the  above  reference  surface  become 
a-,  =  (rQ  +  ?  cos  £2)e0 

a2  =  (c  cos  K2)  t2  +  ^  sin  h\l  ®3 
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and  the  unit  outward  normal  is 


f,  = - * -  [(?  sin  Co)  ,  e„  ■  (?  cos  C0)  9  e,] 

3  +  ?;2  2  -2  r  2  ,z  d 


(3.4) 


where 


er  =  cos  C^  e-|  +  sin  C1  e2 


e0  =  -  sin  C1  e^  +  cos  Cj  e2 


(3.5) 


Next,  certain  basic  properties  of  the  above  reference  surface  are  summar¬ 


ized. 


The  length  of  a  line  element  is  given  by 


ds2  *  aa8  d5a  d?B 


(3.6) 


where 


aa3  =  aa  ‘  aB 


(3.7) 


Thus, 


all  =  (ro  +  K  cos  ^2) 


a22  -  C*  +  ?2>2 


a12  *  0 


(3.8) 


Therefore,  and  Co  are  the  orthogonal  curvilinear  coordinates. 
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The  second  fundamental  form  of  the  reference  surface  is 
da3  •  df  =  baS  dCa  d5s 

where  the  second  metric  tensor  components  are 

{?  sin  K7)  o 

b„  =  -  (r0  +  c  cos  S2)  ^  _ 

J  ^ 

b22  =  -  (C2  +  2  S*2  -  «  2)/^7 
b12  "  0 

Thus,  the  reference  coordinate  systems  are  those  of  the 
curvature  coordinates.  The  radii  of  curvatures  are 

i  .  Vl 

R,  a„ 

1  .  b22 

^2  a22 

The  Gauss-Codazzi  relations  take  the  form 

baB,y  “  bay,8 
or 

d  dr 

f<2  cos  a)  = 

where 

sin  w  =  r/R-| 


(3.9) 


(3.10) 

lines  of 


(3.11) 


(3.12) 


(3.13) 


(3.14) 
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The  non-zero  Crystoffel  symbols,  defined  by 


_  v  —  .  — 

aa,B  =  aB  av  +  baB  a3 


(3.15) 


are  found  to  be 


ril  =  (ro  +  ?  cos  52)(c  cos  h\z^^2  +  ^ 

(?  cos  Z2)  2 

p  1  _  pi  _  _ t-  )  £ 

12  21  “  rQ  +  s  cos  ^ 


(C  cos  S2)  2  ^  cos  ?2)  22  +  ^  sin  h)  22  sin  ^  2 

r22 - 1 - 1 +?r - 1 - *"  <3-16) 


Covariant  differentiation  on  the  surface  is  defined  in  the  usual 


manner, 


3A  .  , 

a  rA  .X 

3?0  ‘  aB  M 

sa"  +  r»  A* 

3Ce  ‘ex  H 


a  -  aB  A  p^  a 

HaB,Y  "  3£y  '  ay  *XB  “  BA  aX 


(3.17) 


where 


Aa  =  aaB  A, 


(3.18) 


The  position  vector  of  a  particle  not  on  the  reference  surface  is 


given  by 
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r  =  r  (Cj ,  K2)  er  +  £3  a3 


(3.19) 


where  £3  is  the  thickness 
The  covariant  base 


coordinate. 

vectors  of  the  shell  space  can  now  be  written  as 


8a. 

9a  =  aa  +  ^3  W 


'a 


g3  =  a3 


(3.20) 


The  metric  tensor  of  the  shell  space  becomes 


a  „  +  2  S,  b  0  +  bj  b  , 

aB  ’3  aB  3  B  aX 


ga3  =  0 
g33  =  1 


(3.21) 


In  the  above  lines  of  curvature  coordinate  system. 


gi  =  ai  n  +  R-) 


g2  =  a2  ^  + 


(3.22) 


Thus,  if  the  thickness  is  small  compared  to  the  radii  of  curvatures, 
then  it  may  be  assumed  that 


gaB  =  aaB 
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det(ga6)  =  det(aaS)  =  a 


(3.23) 


This  approximation  means  that  the  Christoffel  symbols  for  the 
space  are  equal  to  their  value  on  the  reference  surface. 


B.  Field  Equations 

Consider  now  that  the  shell  is  deformed  due  to  external  loads  so 
that  each  point  undergoes  a  displacement 


u  =  u  ?  +  wa3  +  C  i  a0 
a  3  a 


(3.24) 


where  ua  and  w  are  the  reference  surface  displacements  and  wa  denotes  the 
rotation  of  the  normal  to  the  reference  surface. 

Then,  the  Green's  strain  tensor,  according  to  the  approximation 
embodied  in  Equations  (2. 26), takes  the  form 


eaB 

1 

“  2 

<Ua,B 

+  U  „ 

B,a 

kaB 

1 

“  2 

+  w0  ) 
8, a' 

1 

"  2 

K* 

w  -  b 
,a 

The 

stress  resultants 

na* 

-  / 

aB 

T  I 

aB 


.v 


'.a 

1  /  lV 


(3.25) 


a8  t  aB 

m  -  J  t  p  d?, 

h  J 


,a  =  £  TaB 


d£, 


(3.26) 
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Using  Equations  (3.25)  and  (3.26),  the  principle  of  virtual  dis¬ 
placements  yields  the  following  equilibrium  equations: 

naB  +  b“  mvB  +  b“  rv  +  f“  =  0 

,P  v  , p  V 


maP  _  a  _  n 
m,B  r  '  0 


ra  -  naB  b  .  -  maB  bv  b  .  ♦  (nvB  w  ,)  +  p  -  0 

,a  aP  a  vP  ,6  ,a 


(3.27) 


The  corresponding  edge  tractions  take  the  forms: 
nB  =  (naB  +  bB 


a  _  aP 
m  -  m  v, 


r  naB  w_a 


(3.28) 


Finally,  within  the  framework  of  the  initial  stress  formulation  out¬ 
lined  in  Section  II, the  linearized  field  equations  analogous  to  Equations 
(2.26),  (2.28),  and  (2.29)  are 

ect3  ’  \  (ua,B  +  u6,a  +  2  ba8  w) 

kaB  =  \  ("a,B  *  “b,cc>  +  ba  bv8  W  +  ?  (ba  “v,B  +  bB  Uv,a> 


Y  +  w  -  b  uv) 

2  v  a  ,a  av  ' 


(3.29) 
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*  ">*  ♦  "S ^  ’  0 


aB  „a  _  n 

,6  '  '  ° 


“  -  nae  b  .„  -  n.aB  b„0  *  (n"“B  w  „>  +  p  -  0 


,a 


aB 


a  vB 


,  B  i« 


(3.30) 


and 


B  /  aB  .  kB  mavx 
n  =  (n  +  bv  m  )vv 


a  aB 
m  =  m  v. 


a  ,  „0aB 
r  =  r  va  +  n  Wj0i 


(3.31) 


where  n°ag  are  the  initial  membrane  stress  resultants. 

Note  at  this  point  that  if  one  is  to  pursue  the  hybrid  stress  fi¬ 
nite  element  formulation,  then  one  must  satisfy  the  equilibrium  equations 
(3.30)  identically,  which  may  introduce  a  certain  degree  of  numerical  com¬ 
plexity  when  one  deals  with  "deep"  shell  elements.  Certain  simplifications 
may  arise,  however,  when  the  shell  element  is  sufficiently  shallow. 

First  of  all,  recall  that,  for  surfaces,  the  only  non-zero  components 
of  the  Ri emann-Chri stof fel  tensor  RagvX  are  equal  to  ±  R1212-  Moreover  fio] 


R1 21 2  =  ka 


(3.32) 


where  k  is  the  gaussian  curvature  of  the  surface, 


30 


From  the  Gauss-Codazzi  relations, 

R1212  =  b  ;  b  =  det^bct3^  (3.33) 

so 

k  =  b/a  (3.34) 

Consequently,  whenever  k  is  small  enough  or  zero  (developable  sur¬ 
faces),  the  Riemann-Christoffel  tensor  can  be  taken  equal  to  zero, and  then 
the  order  of  covariant  differentiation  on  the  surface  may  be  interchanged. 
In  this  case,  a  stress  function  may  be  employed  to  satisfy  the  homogeneous 
equilibrium  equations.  Otherwise,  recourse  should  be  made  to  the  static- 
geometric  analogy  0  tl- 

For  flat  elements  or  for  developable  surfaces,  one  recovers  the 
plate  equations  of  Section  II. 

C.  Inertial  Loads  Due  to  Rotation 

The  position  vector  of  a  generic  point  in  the  deformed  state  is 

R  =  (r  +  a3)  +  (ua  aa  +  wa3  +  £3  wa  aa)  (3.35) 

If  the  shell  rotates  with  a  constant  angular  velocity, 
fi  =  ft  ?3  (3.36) 

Then,  the  absolute  velocity  of  a  generic  point  is 

v  =  n  x  R  +  (3.37) 
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and  the  acceleration  is 


a  =  f!x(flxR)+2S2x^  +  (3.38) 

For  the  static  analysis, 

#-0  (3.39) 

ot 

so  that 

a  =  Jl2(e3  x  [e3  x  (r  +  S3  a3)  +  e3  x  (u01  aa  +  wa3  +  £  wa  aa)])  (3.40) 

3 

Since  |r|  »  (ua;  w;  53),  the  final  acceleration  is 

a"  =  ft2[e3  x  (e3  x  r)]  (3.41) 

Using  Equation  (3.1),  the  above  acceleration  becomes 

a"  =  -  ft2  xQ  ea  (a  =  1,  2)  (3.42) 

and  the  corresponding  inertia  force  is 

f  =  -  n2  x  e  /  y(S*)d5,  (3.43) 

a  a  (h)  33 

where  y  denotes  tne  mass  density  distribution  across  the  shell  thickness. 
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SECTION  IV 

MATERIAL  CHARACTERIZATION 


A.  One-Ply  Systems 

A  unidirectional  ply  of  cord  and  rubber  is  shown  in  Figure  2.  The 
cord  direction,  denoted  by  1,  and  the  direction  perpendicular  to  the  cord, 
denoted  by  2,  form  a  set  of  directions  referred  to  as  the  principal  di¬ 
rections  of  the  ply.  The  constitutive  law  of  this  cord-rubber  ply  com¬ 
posite  may  be  characterized  by  five  elastic  constants: 

E1  =  longitudinal  Young's  modulus, 

E2  =  transverse  Young' s  modulus, 
v12  =  major  Poisson  ratio, 

\>2i  =  secondary  Poisson  ratio,  and 

G  =  in-plane  shear  modulus. 

Thus,  the  single -layered  composite  has  the  constitutive  law 


e 


12 


(4.1) 


If 


E1  V21  =  E2  v12 


(4.2) 
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Figure  2.  Unidirectional  Ply 


Figure  4.  Simplified  Model 
for  Contraction 


then  the  cord-ply  system  is  referred  to  as  a  structurally  orthotropic 
composite. 

The  rubber  matrix  and  the  reinforcing  cord  are  assumed  to  be 
isotropic.  Thus,  for  the  rubber. 


R  -  1  tnP  »  J*\ 
1  "  Er  ^  1  ‘  R  1’ 


e2  =  (a2  '  VR  al> 


_R 

R 

e12  "  2Gr 


(4.3) 


and,  for  the  cord. 


C  1  ,  C  C> 

=  <01  VC  °2 

C  1  /  C  Cx 

e2  “  Ec  *a2  "  VC  °2) 


C 

_c  _  !]2 
e12  2GC 


(4.4) 


In  the  next  section,  the  elastic  constants  of  the  composite  will  be 
determined  in  terms  of  the  constituents.  It  will  be  assumed  that  the  rub¬ 
ber  matrix  and  the  reinforcing  cords  are  firmly  bonded  together. 
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B.  Tensile  Test 

To  determine  the  modulus  of  elasticity  in  the  fiber  direction,  con¬ 
sider  a  test  specimen  as  shown  in  Figure  3.  The  state  of  stress  of  this 
composite  is 

o?1  ? 1  0 
$2  ~  0 

a]2  =  0  (4-5) 


Thus, 


e2  =  '  v12  E1 

5,2=0  (4-6) 

Assuming  perfect  bond  between  the  rubber  matrix  and  reinforcing 

cords. 


The  condition  of  eaui librium  in  the  1  direction  is 
C  R 

/  a,  da  *  /  da  +  /  a,  da 
A  1  Ac  1  Ar 


where 
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(4.8) 


A  =  total  cross-sectional  area* 

=  total  cord  area,  and 

Ar  =  cross-sectional  area  of  the  rubber  matrix. 
If  one  deals  with  mean  values, 


0|  =  r  o|  +  (1  -  r)c^ 


(4.9) 


where 


The  state  of  stress  in  the  constituents  is 


a!f  t  0 

=  0 

a*  t  0 


a)  =  0 


(4.10) 


<4  =  ° 


(4.11) 


Thus, 


e2  =  "  VR  el 


r 


Mi 


-A 


m 


1 


R  _  n 
e12  '  0 


VC  el 


(4.12) 


Using  Equations  (4.6)  and  (4.12)  in  (4.9),  one  obtains 


E1  =  (1  -  r)ER  +  r  Ec 


(4.13) 


To  determine  the  Poisson  ratio  v12.  consider  the  simplified  model 
shown  by  Figure  4.  Based  on  this  model,  the  total  contraction  is  the 


of  the  contraction  of  the  constituents. 


e2  =  r  e2  +  (1  "  r^e2 


(4.14) 


Putting  Equations  (4.6)  and  (4.12)  in  (4.14) 


v12  =  r  vc  +  (1  -  r)vR 


(4.15) 


C.  Strip  Test 

The  strip  test  in  the  1  direction  is  a  pure  homogeneous  strain,  as 


in  Figure  5, 


e2  =  0 


(4.16) 
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SLIDING  HOOKS 


SLIDING  HOOKS 


Figure  5.  Strip  Test 


Sheet  Test 


Assuming  perfect  bonding, 


o,  =  0 


r)o!;  +  r  a. 


(4.17) 


(4.18) 


0 


(1  -  r)e 2  +  r  4 


(4.19) 


For  the  composite,  one  obtains 


°2  =  E^  V12  al 
el  "  e}  (1  "  V12  V21^ 


(4.20) 


Using  Equation  (4.18)  and  the  averaged  equilibrium  condition  in  the  2  di 
rection. 


°2  =  °2  =  °2 

the  constitutive  relations  for  the  constituents  yield 

el  =  h  (al  ‘  VR 

R  1  /  R\ 

4  ■  rK  (“2  -  al 

E1  -  E-  -  VC  °2> 

e2  =  l~  (°2  '  VC 

u 


(4.21) 


(4.22) 


(4.23) 
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From  Equations  (4.20),  (4.22),  and  (4.23), 


a?  =  ER  £1  +  VR  v12  al 


al  =  EC  el  +  VC  E^  v12  al 


(4.24) 


The  substitution  of  Equations  (4.24)  into  (4.1)  yields 


E2  v12  =  E1  V21 


(4.25) 


Using  Equations  (4.20)  and  (4.24)  in  (4.22)  and  (4.23),  one  obtains 


R  1  -  VR 

e2  =  ~~T^~  V21  al  '  VR  el 


c  1  ■  VC 

e2  =  E„  V21  al  ‘  VC  el 


(4.26) 


Putting  Equation  (4.26)  in  (4.19),  and  employing  (4.15),  yields 


v 


°1  v?1  (1  -  r)(l  -  v£)Fc  +  r(l  -  v£)Ef 


12 


ER  EC 


(4.27) 


Comparing  the  above  expression  with  Equation  (4.20),  one  obtains 


‘I 


12 


er  ec 


1  -  v12  v2;  ‘  V21  (1  -  <•)(!  -  v|)Ec  +  r(l  -  v^)Er 


(4.28) 


With  Equation  (4.25),  the  above  expression  yields 


er  ec 


1  -  V,0  v 


12  21 


(1  -  r) (1  -  ^)EC  +  r(  I  -  v‘)tf 


(4.29) 


Dr 


er  ec  ei 


r) (1  -  v^)Ec  +  r(1  -  v£)tRJ  +  ER  Ec  v(2 


(4.30) 


Eliminating  E]  and  u,2  from  Equation  (4.30),  one  obtains  the 
secondary  Young's  modulus  in  terms  of  the  constituent  properties, 


Er  Ec  [(1  -  r)ER  +  r  Ec] _ 

E2  =  Er  Ec  +  r(l  -  r)L(ER  -  £<0“  -  (ER  v'c  -  Ec  vRPJ 


(4.31) 


D.  Sheet  Test 

The  sheet  test  is  a  pure  homogeneous  strain  with  (Fig.  6) 

(4.32) 

el  =  e2  =  e 

Again,  employ  the  perfect  bond  condition  and  the  equilibrium  in  the  1  and 
2  directions, 


e  =  e2  0  '  r)  +  r  e2 


o1  =  (1  -  n)a^  +  r 
a2  =  a2  =  a2 


(4.33) 
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For  the  composite,  one  obtains 


P  _  !l  vjin 

e  "  e  "  E  °2 

bi  H  c 


VJ2n 

e  ~  E  ”  E  al 

t2 


(4.34) 


or 


a,  = 


1  1  -  v)2  v21 


(4.35) 


For  the  constituents, 


1-  (a*  -  vR  a2) 


^  (a2  '  VR  al} 


(4.36) 


and 


el  =  Ec  (al  "  VC  °2^ 


e2  =  Ec  (a2  "  VC  al^ 


(4.37) 


Following  the  steps  of  the  previous  section,  one  obtains 


E2  v12  =  E1  V21 


(4.38) 


and 
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(4.39) 


E.  Pure  Shear 

Consider  the  test  illustrated  in  Figure  7.  At  the  center  of  the 
specimen,  the  state  of  stress  may  be  characterized  by 


From  static  equilibrium, 

R  -  rrC  -  a 
a12  '  12  "  12 

Thus,  the  constitutive  relations  read, 

°12 
e12  "  G 

R  _  °J2 
el2  ‘  2Gr 

C  _  012 
el2  “  2EJT 


(4.42) 


(4.43) 


•tmn 


Using  Equation  (4-43)  in  (4-41),  one  obtains 

G  GCGR 

(1  -  r)Gc  +  r  Gp 

F.  Summary  of  In-Plane  Test  Results 

The  elastic  moduli  for  the  one-ply  composite  are: 

E,  =  Ec(l  -  r)X  +  r 

v,2  -  r)vR 

r-  _  (1  -  r)X  +  r 

h  '  ER  X  +  r(1  -  r)[T5i  -  l)2  -  (Xvc  -  vR)^] 

E2 

v2i  ■  rj-  vi2 
G  =  GC  (1  -  r>  +  m 

where 

X  =  er/ec 
a)  =  gr/gc 


(4.44) 


(4.45) 


(4.46) 


Note  that  Equation  (4.45)  is  not  an  a  priori  assumption,  but  a  con¬ 
sequence  of  this  modeling  approach. 

These  formulas  should  be  compared  with  the  experimental  results  of 
the  previously  prescribed  tests. 
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Thus,  in  the  tensile  tests, 

(a)  measure  Cp  ,  and  £2,  and 

(b)  calculate  E-|  and  v^. 

In  the  strip  or  sheet  tests, 

(a)  measure  a ^  and  ,  and 
(d)  calculate  E^/(l  -  v21)  =  a]/£i> 

from  which 


v 


al  "  E1  el 


21 


12 


Finally,  using  Equation  (4.45), 


E„  = 


E1  V21 


2  v 


12 


G.  Material  Properties  in  the  Transverse 
Direction 


In  the  longitudinal  direction, 


Tio  =  r  +  (1  "  r)T 


13 


R 

13 


t13  =  2  Gc  e13 


t13  =  2  GR  e13 


t13  s  2  G13  e13 


Thus, 


G1 3  =  r  GC  +  (1  "  r)GR 
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(4.47) 


(4.48) 


(4.49) 


(4.50) 


e23  "  ^  “  r^e23  +  r  e23 

R  C  _  _ 
t23  "  x23  "  t23 

g23  =  T23/2GR 

g23  =  t23/2GC 

e23  =  t23^2G23 


(4.51) 


(4.52) 


H.  Transformation  to  Local  Frame 

In  general,  the  local  coordinate  system  does  not  coincide  with  the 
coordinate  system  used  to  formulate  the  layer  constitutive  relations. 

At  each  point  on  the  reference  surface,  the  cord  orientation  may 
be  defined  by  the  vector 


c  =  cos  8 


al  a2 

—  +  sin  6  — 


(4.53) 


where  the  cord  angle  8  is  specified  by  the  lift  equation  [7] 
cos  8  *  (r  +  4  cos  S2)cos  y/rb 

rb  =  bead  radius 

Y  =  Green  angle  (4.54) 
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In  the  principal  frame  of  reference  (3  =  0),  the  constitutive  rela¬ 


tions  have  the  form 
P*  -  !l  V21 

E1  "  F  F 

tl  L2 

*  _  °Z  v12 

2  "  E2  '  E1 

Y12  =  CT1 2/Gl 2 

y13  =  a13/G13 

★  * 

y23  =  a23/G23 


(4.55) 


The  new  elastic  constants  in  the  lines  of  curvature  coordinate 
system  are,  [12] 

ell  =  C11  T11  +  c12  t22  +  c13  t12 

e22  =  c12  T11  +  c22  t22  +  c23  T12 

y12  =  C1 3  T11  +  c23  t22  +  C33  T12 

Y13  =  S11  t13  +  s12  t23 

y23  =  s12  T13  +  S22  t23  (4.56) 


where 


I .  In-Plane  Cord  Angle  Change 

Consider  a  filament  emanating  from  a  material  point  P  in  the  un¬ 
deformed  configuration,  characterized  Ly  the  vector  dTT.  After  deformation, 
the  same  filament  will  be  characterized  by  the  vector  dF,  Figure  8.  If 
F  is  the  deformation  gradient,  then,  [13] 

dF  =  F  dR  (4.58) 

Let  the  displacement  vector  be  lT,  such  that 

u  =  F  -  W  (4.59) 

Thus, 

du  =  dr  -  dR  =  (F  -  I)dlf 

•V  "V 

where  I  is  the  unit  map. 

The  stretch  is  defined  by 

X  =  ^ 
dS 

where 

dS  *  |dR| 

ds  =  | dr |  (4.62) 


(4.60) 


(4.61) 


Using  Equation  (4.58),  one  obtains  the  following  expression  for  the  stretch. 


Xn  =  Fff 

X2  =  N  •  FT  FIT 

«■«* 


(4.63) 


mm&m 
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Figure  8.  Angle  Change  Between  Two  Fibers 


where 


n  = 


dr 

ds 


N  = 


dR 

dS 


Define  the  displacement  gradient  as 


(4.64) 


du  =  A  dR” 


(4.65) 


Comparing  the  above  definition  with  Equation  (4.60),  one  obtains 

A  =  F  -  I  (4.66) 

The  classical  strain  tensor  is  defined  by 

E  ■  \  <FT  F  -  I)  (4.67) 

or,  using  Equation  (4.66), 

E  =  |  (A  +  AT  t  AT  A)  (4.68) 


In  terms  of  the  above  strain  tensor,  the  stretch  defined  by  Equa¬ 
tion  (4.63)  is 

X2  =  2  N  •  E  N  +  1  (4.69) 

Note  that  if  E  =  0,  then  the  stretch  is  unity. 
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For  deformations  in  the  small  strain  domain,  it  is  convenient  to 


define  the  extension  6, 

6  =  X  -  1  (4-70) 

Using  Equation  (4.70),  one  obtains  the  following  expression  for 
the  extension, 

62  +  2  S  =  2  N  •  E  N  (4.71) 

If  S2«  1 ,  then 

5  -  N  •  E  N  (4-72) 

and  it  is  spoken  of  as  small  extension. 

In  component  form,  Equation  (4.72)  reads 

6  =  E  .  Na  N3  (4-73) 

aB 


where 


E  =  Ea6  (5°  x  a*  (4-74) 

where  a°  and  a”  are  contravariant  and  covariant  base  vectors,  respectively. 

a 

Consider  now  two  fibers  in  the  undeformed  configuration,  defined 


C2  =  C2  aa 


(4.75) 


as  shown  by  Figure  8.  The  shear  of  these  two  fibers  is  defined  by  y, 


Y  =  ©  -  9 


(4.76) 


where 


©  =  the  angle  between  the  two  fibers  in  the  undeformed  configura¬ 
tion,  and 

0=  the  angle  between  the  same  two  fibers  in  the  deformed  config¬ 
uration. 

Using  Equation  (4.63),  the  corresponding  deformed  directions  are 


nl  =  ~  ^1 


n2  =  X2  ~  C2 


Thus, 


cos  0  =  n]  •  n2  =  x^  C2  •  FT  F  (T-, 


Since 


F'  F  =  2  E  +  I 


then  Equation  (4.78)  yields 


cos  0  =  x  ~*x-  [2  T2  •  E  +  cos  @] 


(4.77) 


(4.78) 


(4.79) 


(4.80) 
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where 


(4.81) 


cos  ®  =  C.|  • 


Using  Equation  (4.76), 


cos  8  =  cos  ©  cos  y  +  sin  ®  sin  y 


(4.82) 


and  the  component  form  of  (4.80)  becomes 


cos  ©  cos  y  +  sin  ©  sin  y  =  [2  ^  C2  +  cos  ®3  (4-83) 


Note  that  the  usual  (©  =  tt/2 )  shear  is 


1  ro  r  r8- 


sin  Y  =  t— y 
A1  a2 


p  EaB  <1  CP 


(4.84) 


For  small  extensions,  Equation  (4.83)  yields 


cos 


©  cos  Y  +  sin  ©  sin  y  =  2  E^  C1  C2  +  cos  © 


(4  .85) 


Knowing  c“,  C^,  ©,  and  the  strain  tensor  Ea3,  Equation  (4.85) 
is  a  transcendental  equation  for  the  shear  y.  For  small  shear. 


cos  y  -  1 
sin  y  -  y 
one  obtains 


(4.86) 


Y  = 


2  ci  j 

sin  © 


(4.87) 
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Knowing  the  shear  y,  the  angle  between  the  two  fibers  after  deformation  is 


©  =  ®  - 


2  E 


"a  8 


Li  4 


sin  ® 


(4.88) 


The  above  cord  angle  change  is  not  yet  incorporated  into  the  com¬ 


puter  code. 


SECTION  V 

NUMERICAL  FORMULATIONS 

A.  Discretization 

Flat  triangular  elements  are  used  with  nodes  being  located  on  the 
reference  surface.  For  the  rotationally  symmetric  problem,  the  global 
coordinate  system  is  that  of  the  lines  of  curvature  coordinates  of  the 
reference  configuration. 

Each  of  these  triangles  is  defined  by  its  vertex  position  vectors 
r*|,  r^,  and  ry 

Using  the  vertex  position  vectors,  the  following  local  base  vectors 
are  generated, 

?i  = 

g2  =  g3  x  g, 

g3  =  (?3  X  ?1)/|r3  X  Ci  I  (5,1) 

where 

*1  =  ?3  '  ?2 

7  =  f  -  r,  (5-2) 

r2  1 

The  origin  of  this  local  coordinate  system  is  taken  to  be  the  cen¬ 
troid  of  the  triangle.  The  cartesian  components  of  the  local  base  vectors 


are  defined  by 


g.  =  g. j  e^  (i,  j  =  1,  2,  3)  (5.3) 

The  area  integrations  are  carried  out  by  gaussian  cubatures 
over  T2  simplexes  [14].  The  procedure  is  now  outlined  for  an  integral  of 
the  form 


l  -  I  F(n1 ,  n2)dn1  dn2 


(5.4) 


Let  the  in-plane  vertex  coordinates  be  denoted  by  v^  and  va2  in  a  bary- 
centric  local  coordinate  system.  In  order  to  use  gaussian  cubature  formu¬ 
lae  for  the  numerical  integration  of  Equation  (5.4) ,  each  triangle  is 
mapped  onto  a  standard  triangle  (0,0;  1,0;  0,1)  by  the  following  simplex 
transformation: 

Z1<V  n2>  '  al  a1  +  a2  n2  +  a3 

z2(n-|,  n2)  =  b1  n1  +  b2  n2  +  b3  (5.5) 


where  the  constants  a^  and  b.  are  calculated  from  the  transformation  con¬ 


straints: 


and 


z2(vn*  V12^  =  0 
z2 (v2i »  ^22 ^  ~  ^ 

Z3^V31 ’  V32^  =  1 

Then,  using  the  inverse  of  Equation  (5.5), 

z2)  =  [ (z.|  -  a3)b2  -  (z2  -  b3)b1]/(a1  b2  -  a2 

n2(zr  z2)  =  [(z2  -  b3)a1  -  (z-|  -  a3)a2]/ (a-j  b2  -  a2 

the  integral  defined  by  Equation  (5.4)  becomes 
1  1-z 

I  =  |T|  /  /  F[n] (z] ,  z2);  n2(zr  z2)]dz]  dz2 

0  0 


where  the  Jacobian  T  is  defined  by 

3n-|  3n1 

3 (V  ’l?) 

T  =  3T^T^T  = det 


3z.|  3z2 


which  is  calculated  from  Equation  (5.8). 

Next,  the  limits  of  integration  in  Equation  (5.9)  are 
[-1,  1]  by 


z]  =  (1  -  0/2 


(5.7) 


})  (5.8) 


(5.9) 


(5.10; 


changed  to 
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i  =  UL 
1  8 


1 


J 


-1 


f(Z,  n)(l  +  OdC  dn 


(5.12) 


where 


f(Z,  n)  =  FCn^z^O.  z2(S,  n2Cz1(0,  z2(5,  n)]}  (5.13) 


Equation  (5.12)  is  now  in  the  desired  form  to  use  tabulated  gaussian 
weights  and  nodes.  For 


It!  G  G 

I  =  V  E  S  C..  f(E.,  n,) 

8  i=l  j=l  1J  1  J 


(5.14) 


where  the  weights  Ai  and  B,  and  the  nodes  £.  and  are  tabulated  in 
[14]  according  to  the  following  schedule: 

1  G 

/  h ( x ) ( 1  +  x)dx  =  E  A.  h(x.) 

-1  i=1  1  1 

1  G 

/  h(x)dx  =  E  B.  h(x. ) 

-1  i=l 

An  alternative  form  of  Equation  (5.15)  is 


(5.15) 

l 
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(5.16) 


,.JIL  2 


Cfl  f<PSl> 


H  =  G  x  G 

where  is  the  gaussian  node  defined  by  the  pair  (£fi,  n^). 

Product  integrations  are  performed  as  follows: 

[I]  =  /  /  CF(n1 ,  n2)]T  [D]  [G^,  n2)]dn1  dn2  (5.17) 

where  the  compliance  matrix  is  assumed  to  be  constant  over  the  element. 
Thus,  one  way  to  perform  the  designated  integration  is  to  carry  out  the 
multiplication  to  obtain 


[I]  ■  J  /  [E(n-, ,  n2)]dn-|  dn2  (5.18) 

Anin2 

where 

E(n-| ,  n2)  =  FT  (nj,  n2)D  G(n(>  (5.19) 

Each  term  in  the  E  matrix  is  a  simple  polynomial  whose  integral  may 
easily  be  tabulated  using  the  above  procedure. 

However,  at  this  development  stage,  a  more  flexible  approach  is 
adopted  which  performs  integration  on  matrix  products  in  general.  Thus, 


u 

[I]  =  M  ^  Clj  Cf(ri.  rij)]  [D]  [g(s1(  n.)] 


(5.20) 


i,j=l 
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B.  The  Hybrid  Stress  Finite  Element 
Formulation 

The  hybrid  stress  functional  for  the  problem  under  consideration  has 
the  form,  from  Equation  (2.43), 

N  l  T 

*  =  1  2  //{/{^(v  '2;  x3^  CD(x3>DCcr(n-| »  n2;  x3)] 

n=1  rn 

+  2[x(n1,  ty,  x3)]T  [D(x3)][t(ti1  ,  ry  x3)])dx3>dn1  dn2 

-  j  [N(s)  •  u(s)  +  M(s)  •  ^(s)]ds  (5.21) 

RN 

where 

D  =  direct  stress  compliance, 

S  =  shear  stress  compliance, 
a  =  column  (cy,  a 22>  a^)  =  direct  stress, 
t  =  column  (a^3,  a23)  =  shear  stress, 

N(s )  =  (Ny ,  Nt,  A)  =  boundary  moment  resultant  vector, 

M(s)  =  (M^ ,  M  )  =  boundary  moment  resultant  vector, 
tT(s )  =  boundary  displacement  vector,  and 
<j>(s)  =  boundary  rotation  vector. 

The  stress  field  for  each  layer  is  chosen  such  that  the  layer  equil¬ 
ibrium,  interlayer  equilibrium,  and  overall  resultant  equilibrium  for  an 
element  are  satisfied.  However,  stress  continuity  along  the  interelement 
boundary  need  not  be  satisfied.  Accordingly,  the  solution  of  the  equilib¬ 
rium  equations  is  obtained  as  the  sum  of  the  homogeneous  and  particular 
solution,  such  that  the  homogeneous  state  of  stress  satisfies 
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croo  +  0o,  =  O 
a8,8  u3,3 


a  0  =0 

a3,a 


(5.22) 


throughout  the  thickness  of  the  shell.  Thus,  for  each  layer  i, 

a1OQ  +  a\,  =  0 
a8,8  a3,3 

aa3,a  =  ®  (i=l,2...L)  (5.23) 

where  L  is  the  number  of  layers.  In  order  to  satisfy  interlayer  equilibri 
one  must  have 

aa3^T  n2’  hi  ^  =  aa3  ^nl  ’  n2’  hi  ^  ^5'24^ 

where  h^  is  the  thickness  coordinate  of  the  interface. 

The  particular  stress  field  is  obtained  in  terms  of  the  stress  re¬ 
sultants  which  satisfy  the  incremental  equilibrium  equations 


NkR  +  Fk  =  0 

a3,8  a 


MaB.e  '  i  -  0 


Qk  +  (Nk_1  wk  )  +  Fk 

ya,a  '  ay  ,y  ,a  3 


=  0 


(6.25) 


k  k 

for  the  k-thstep  in  an  initial  stress  formulation,  where  Fq  and  are 

the  k-th  load  increments,  and  w  is  the  current  normal  displacement.  In 

k 

the  case  of  the  rotating  tire,  the  external  loads  F^  are  replaced  by  the 
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1/ 

associated  inertial  forces.  After  calculating  the  stress  resultants  N 

a3 

k  k 

Mag,  and  Qa  in  terms  of  the  incremental  loads,  the  particular  stress  field 
is  obtained  in  the  usual  fashion, 


4 2  h  I’  • 4  (fr)’  *4  <5-26> 

where  h  is  the  shell  thickness. 

Thus,  the  total  stress  field  for  each  layer  will  have  the  form 


^  ~  ag(T1i »  ^2*  x3^  °a(V  q2  ’  x3^a  aq(rii  *  ^2’  x3^ 

T  =  W  n2’  X3^6  +  Ta^r  n2’  x3^a  +  aq^nr  n2;  x3^q  (5*27) 


where 

a6(nr  n2;  x3)3  =  homogeneous  direct  stress  field  , 

aot(nl  ’  n2’  x3^a  =  Part1’cular  direct  stress  field  due  to  the  exter- 
'  c  J  nal  incremental  loads,  and 

0q{nr  n2'»  xjq  =  particular  direct  stress  field  due  to  the  presence 
of  normal  displacements  in  the  third  resultant 
equilibrium. 

with  identical  definition  for  the  shear  stress  distribution  x.  Further¬ 


more, 

3  =  column  (3-|,  32,...32g)  =  undetermined  stress  coefficients, 

a  =  column  (a^ ,  ot2, - a,g)  =  external  force  vector  obtained  by 

linear  interpolation  of  the  nodal  values, 
and 
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q  =  column  (q,,  q«,...q1(.)  =  nodal  displacement  vector  for  an 
L  5  element. 

The  vectors  a  and  q  are  common  to  all  layers,  while  the  vector  3 
for  each  layer  is  such  that  the  interlayer  equilibrium  (5.24)  is  satisfied. 

For  each  element,  the  edge  resultants  may  be  calculated  from  the 
stress  field  (5.27)  to  yield 

l?(s)  =  [N®(s)B  +  N^(s)a]v  +  [NB(s)  +  N“(s)]t  +  [QB(s)S 
+  Qa(s)a  +  Qq(s)q]n 

M(s)  =  [MB(s)3  +  M^(s )a  +  Mq(s)q]v  +  [MB(s)g  +  M^(s)a 

+  Mq(s)q]t  (5.28) 

where  v  and  t  are  the  outward  normal  and  tangent  vectors  along  the  edge 
of  an  element  and  n  is  the  outward  normal  to  the  plane  of  the  element. 

The  displacement  and  rotation  vectors  along  the  boundary  of  an 
element  are 

u(s)  =  (Ut(s)q)t  +  (Uv(s)q)v  +  (W(s)q)n 

♦  (s)  =  (4>t(s)q)t  +  (<J)v(s)q)v  (5.29) 

where  the  vectors  Ut(s),  Uv(s ) ,  W(s),  d>t (s )  and  4>v(s)  imply  interpolation 
along  the  boundary  of  an  element  in  terms  of  the  nodal  displacements  q-j , 

q2 ’ ‘ ' ql 5‘ 

Using  Equations  (5.27),  (5.28),  and  (5.29)  in  (5.1),  one  obtains 
the  following  hybrid  stress  functional  suitable  for  numerical  computation. 
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N  1  T  T  T 

tt(6,  q)  =  ^  {1  [3  CgB  3  +  2  31  PB  +  2  31  AC6q  q 


n=l 


Bq 


+  2  aP^  q  +  qT  ACqq  q]  -  [f5T  «Bq  q  +  zj  q  +  qT  AWqq  q])  (5. 


where  the  vector  3  refers  to  all  the  stress  coefficients  for  an  element, 
and  where 


L 


C33  =  i^1  ^c3S^i 


p3  =  1;1  <pB>i 


AC 


Bq  =  i”  (ACBq}i 


L 


AP  =  l  (Ap  ) . 
a  -j-]  a  1 


AC  =  l  (Ac  ) . 

qq  i=1  v  qq'i 


WBq  *  ,!?,  (Vl 


za  =  ”aq  a 


with 


h.. 


CSB  ’  (  J  [o6  +  2  TB  SVdx3  dnl  dn2 

A  n.+] 
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30) 


'EE: 


>  '  M<*»T  UX)<S>  +  [Nt<S”T  Ut(5)  +  [Mt(s)]T  ^<*> 

+  [M®(s)]T  *t(s)  +  [qs(s)]T  w(s)}ds 

“aq  *  ^  {[N“!s)]T  Uv(s)  +  K(S)]T  Ut(s)  ♦  [MJ(s)]T  *v(s) 

+  [^(s)]T  *t(s)  +  [f(s)]T  VJ(s )>ds 

AW  q  -  j  {[mJ(s)]T  »v(s)  +  [Mq(s))T  *t(s)  +  [Qq(s)]T  W(s)}ds  (5.31) 

For  clarity  in  presentation,  the  following  nomenclature  is  intro¬ 
duced: 
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Cgg  =  element  flexibility  matrix, 

Pg  =  element  flexibility  vector, 

ACgq  =  homogeneous  incremental  element  flexibility  matrix, 

APa  =  incremental  element  flexibility  vector, 

ACqq  =  particular  incremental  element  flexibility  matrix, 

WBq  =  hybr"'c*  eTement  load  matrix, 

Za  =  hybrid  element  load  vector,  and 
AWqq  =  incremental  hybrid  element  load  matrix. 

The  symbol  M  in  Equations  (5.31)  implies  merging  procedures  accord¬ 
ing  to  the  interlayer  equil ibrium  condition  (5.24). 

Since  the  stress  coefficients  3  are  independent  for  each  element, 
the  variation  of  tt  with  respect  to  B  yields  for  each  element, 

C$8  8  +  P8  +  %  "  *  %  *  0  <5'32> 
from  which 

8  *  -  C8B  [pe +  (sceq  -  Vq]  (5-33) 

Next,  substituting  Equation  (5.33)  into  (5.30),  and  taking  the  variation 
with  respect  to  q,  one  obtains 
N 

Z  (kr  -  P)  =  0  (5.34) 

n=l  n 

where  k  and  p  are  the  element  stiffness  and  load  matrices,  defined  by 
k  =  k1  +  k2 

P  =  P1  +  P2  (5.35) 
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where 


k2  '  ACeq  CB8  AC8q  +  *ACSq  C6B  W8q*  +  ^CBq  CB8  W8q' 

+  ACqq  +  ?  (AWqq  +  AWJq> 

pi  =  za  '  “jq  cbs  ps 

P2  -  -  APa  *  AcT6q  CS8  Pe 

where 

k1  =  linear  hybrid  stress  element  stiffness  matrix, 

=  incremental  hybrid  stress  element  stiffness  matrix, 
p-j  =  linear  hybrid  stress  element  load  vector,  and 
P2  =  incremental  hybrid  stress  element  load  vector. 

From  this  point,  the  analysis  follows  the  usual  steps  of 
placement  method.  Thus, 

Kr  =  g  . . . 

where 

N 

K  =  Z  (k)  =  structure  stiffness  matrix 
n=l  n 

N 

g  =  Z  (p)  =  structure  load  vector 
n=l  n 

r  =  unrestrained  generalized  displacements 


(5.36) 


the  dis- 


(5.37) 


(5.38; 
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Finally,  Equation  (5.37)  is  solved  for  the  unknown  displacements  r 
which  are  used  to  update  the  geometrical  configuration  and  to  calculate 
the  initial  stress  resultants  exhibited  in  the  resultant  equilibrium  equa 
tions. 

B.  Displacement  Formulation  with  an  Assumed 

FTicewise  Parabolic  Shear  Stress  Distribution 

The  linearized  field  equations  of  Section  II  are  first  summarized 

here  for  clarity  in  presentation. 

(a)  strain-displacement  equations: 
ea3  =  1  (Ua,B  +  U3,J 
^a3  7  ^Wa,(3  + 


Y  =  oj  +  w 

'a  a  ,a 

(5.39) 

(b)  equilibrium  equations: 

naB,B  +  fa  =  0 

maB,B  -  ra  *  0 

ra,a  +  (naS  w,c>,B  +  P  =  0 

(5.40) 

(c)  boundary  conditions: 


on  S^,  and 


"a  "  V  vs 


ma  =  maB  VB 


r  *  ra  va  +  "aB  “,a  VB 


on  S-|. 


(d)  strain-stress  relations: 


dD  /  \ 

ectB  =  9n^  na£’  ma6*  V 


kaB  =  3^  (naB’  "W  ra> 


9r  ^naB*  maB*  ra^ 


The  Reissner  functional  for  the  above  field  takes  the  form 


"R  ’  l  {-  B<’W  ma3*  rc)  +  -f  w,cx  “,B  ‘  fo  ucc  '  pw 

+  “T  (ua,B  +  u6»a)  +  "T  K.,6  +  “B.a>  +  ra(“a  +  w,a),dfl 


-  L  <"a  ua  +  "a  "a  +  ™>dS 


I  [na(ua 


+  mfw  -  wj  +  r(w  -  w)]dS 


The  variational  equation  6tt^  =  0  yields 

(a)  equilibrium  equations, 

(b)  displacement  boundary  conditions, 

(c)  traction  boundary  conditions,  and 

(d)  displacement  gradient-resultant  relations. 

This  variational  principle  will  now  be  utilized  to  relate  the  ref- 

k  k 

erence  surface  strains  to  the  layer  stresses  Tag  and  where  k  denotes 
the  layer  number.  First,  in  order  to  cope  with  an  arbitrary  shear  stress 
distribution  across  the  thickness,  a  piecewise  parabolic  shear  stress 
field  is  introduced  in  the  following  form: 

(5-45) 

and  the  direct  stresses  will  be  represented  by 

TaS  =  Nk(z)naS  +  Mk(z)l"aS  (5'46) 

where  k  indicates  the  k-th  layer,  and  the  functions  Mk(z),  Nk(z),  and 
1/ 

S  (z)  will  be  determined  later  in  Section  V.C.l. 

For  the  purposes  of  the  present  discussion,  it  is  sufficient  to  note 
that  S  (z)  is  a  second-degree  polynomial  in  the  thickness  (z)  coordinate. 

In  each  layer,  the  local  strain-stress  relations  are  (see  Figure  9 
for  cord  angle  description) 


k  _  _k  k  .  k 


r\ 

11  = 

_  r\  r\ 

C11  T11 

+  r 

22  c13 

x12 

k 

k  k 

k 

k  ,  k 

k 

22  " 

c12  T11 

c22 

t22  +  c23 

x12 
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k  _  _k  k  k  k  +  k  k 

Y12  c13  11  c23  22  33  12 


f13 


k 

y23 


sk  Tk  +  k  k 
S11  13  s12  t23 


k  k  .  k  k 
s12  T1 3  s22  t23 


(5.47) 


where  the  constants  ckj  and  skj  are  given  by  Equations  (4.57)  for  the  layer 
under  consideration. 

In  short. 


ek  =  Ck  xk 
aB  aBY<5  Y<$ 


k  _  nk  k 

Ya  aB  33 


(5.48) 


The  complementary  energy  density  for  the  k-th  layer  is 


pk  _  pk  k  k  pk  k  k 

aBY<5  otB  yS  aB  a3  B3 


(5.49) 


n 


_ 

‘j 


* 


Using  Equations  (5.45)  and  (5.46)  in  (5.49),  one  obtains  the  com¬ 
plementary  energy  density  in  Equation  (5.44)  as  follows: 


B(,W  "W  ra> 


L  k+1  k 

kf,  J  B  <%!'  V  V  z)dz 


(5.50) 


It  follows  that  the  variationally  consistent  reference  surface 
strain-resultant  relations  are  those  of  Equation  (5.43).  The  inversion 
of  these  relations  yields  the  resultant-strain  relations  in  the  form: 


naB  =  Aa3Y6  EY<5  +  BaBY$  kT<5 


mct8  ^a8y5  +  ^a8y<5  ky6 


ra  =  Sa8  Y8 


(5.51) 


At  this  stage,  recourse  is  made  to  the  potential  energy  functional 


V([E(V  "•“0)+¥w,aw.8-faUB-'),']dA 


-  /  (n  u  +  rn  w  +  rw)dS 

L  v  a  a  a  a 


(5.52) 


and  thus,  the  analysis  will  be  based  on  the  displacement  variational  prin 
ciple. 


1.  Piecewise  Parabolic  Shear  Stress  Distribution 

For  each  layer,  the  components  of  the  stress  tensor  are  expressed 
in  terms  of  the  stress  resultants  as  follows: 


k  =  1  d2  ek 
TaB  A  dz* 


1  d  Qk 
B  dz 


-  ek  (z)r 


a 


(5.53) 


where  k  denotes  the  layer  number.  Figure  10,  and  the  constants  A  and  B, 
and  the  functions  6  (z)  are  established  such  that 

L  Ck+1  k 

nce%_£,  l  Wv  n2;  z)dz 

K~  I  Ci. 
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1  z  4(T1]’  n2;  z)dz 

'  ck 

L  Ck+1  k 

r  =  z  !  t*  (nr  n2;  z)dz  (5.54) 

a  k=l  ck  aJ  1  L 

where  ck  and  ck+1  denote  the  k-th  layer  surface  positions  as  measured  from 
the  reference  surface,  Figure  10. 

For  each  layer,  the  usual  three-dimensional  equilibrium  equations 
are  to  be  satisfied, 


k  .  k  ,  *k  n 

Tnn  +  T-  +  f  =  0 

oc6,3  a3,z  a 


4,c. +  <4  ",b>,o.  +  4,z 


(5.55) 


The  overall  equilibrium  requires  that  the  shear  stresses  are  continuous 
across  the  boundaries  of  adjacent  layers, 


Tot3^ck+l ^  =  Ta3  ^Ck+1 ^  =  1 »  2, . . . (L-l )] 


(5.56) 


On  the  inner  and  outer  surfaces,  the  shear  stresses  are  zero, 
Ta3^cl ^  =  0 

Ta3^CL+l ^  =  0  ,57^ 

Equations  (5.54)  require  that 

A  -  Z  r  (z)dz  (5.58) 

k=!  ck 
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- . . 


i-nwmmmmmi 


L  ck+l 


0  =  E  / 
k=1  ck 


(z)dz 


(5.59) 


L  k+1  H2«k  , 

0  =  E  f  z  (z)dz  (5.60) 

k=1  ck 


_  L  •  k+1  j  Jc 

B  =  Z  /  z  4j§_  (z)dz  (5.61) 

k=1  ck 

L  ck+l  k 

B  =  -  E  I  ©k(z)dz  (5.62) 

k=1  ck 

Now,  Equations  (5.56)  and  (5.57)  imply  that  Equation  (5.59)  is 
identically  satisfied.  It  also  follows  that  Equation  (5.61)  is  identical 
to  that  of  Equation  (5.62). 

Equation  (5.60)  essentially  defines  the  reference  surface  location. 
Figure  10. 

If  the  shear  stress  distribution  across  the  shell  thickness  is 

K 

approximated  by  a  piecewise  parabolic  distribution,  then  0  (z)  has  the 
form 

©k(z)  =  sk(z2  -  c£)  +  Sk  (5.63) 

where  sk  and  Sk  are  characteristics  of  the  layer  under  consideration. 
Using  the  above  form, 

A  '  W  2  Sk<Ck+1  '  Ck) 
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B  =  3  sk(ck+l  “  ck^ 


kf1  sk(ck+l  "  ckJ  =  0 


From  Figure  10, 


ck  ’  bk  -  h 


so  that  the  third  of  Equations  (5.64)  takes  the  form 


h  =  2  kf1  5k(bk+l  "  bk)/kf1  sk^bk+l  "  M 


From  Equations  (5.56)  and  (5.57),  one  obtains 


s1  =  0 


Sk  =  Sk-1  +  sk-l^ck  '  ck-l  ^  =  2’  3”*>L) 


The  stress  tensor  components  for  the  layers  are 


k  _  2sk 


2s, 


_  r\._  i  in  _ 

ToS  '  TT  naB  T  Z 


a8 


Ta3  =  -  1  [sk(z!  -  ck>  +  Sk]ra 


Note  that  f*  in  Equation  (5.55)  can  be  shown  to  be 

k  2sk 
fK  =  -ir-  f 

a  A  a 


(5  764) 


(5.65) 


(5.66) 


(5.67) 


(5.68) 


(5.69) 


and 


P  =  T33(CL+1>  '  T33<c1> 


(5.70) 


The  transverse  normal  stress  from  the  second  equations  of  (5.55) 
and  (5.68)  is  found  to  be 


4<z> =  f  (sk[i (z’  -  cP  -  ck(z  -  ck)]  +  sk(z  •  ck> 


•k'  k 


+  Rk}ra,c. 


(5.71) 


in  the  absence  of  the  initial  stress  term,  where  is  chosen  such  that 
the  normal  transverse  stress  is  continuous  across  the  boundaries  of  adja¬ 
cent  layers. 

Using  Equations  (5.68),  a  variationally  consistent  set  of  reference 
strain-resultant  relations  are  now  derived  according  to  the  outline  of 
Section  V.C. 

For  compactness,  the  following  notation  is  introduced: 


nl  =  rV  el  =  ell 

u2  =  n22  e2  =  e22 


81 


I-1 


n3  =  n12 


n5  =  m22 


n6  =  m12 


n?  =  r1 


n8  =  r2 


e3  =  2  £12 


-H  “i-i— 

e5  =  k22 
e6  =  2  k12 


e7  =  Y12 
e8  =  Y13 


(5.72) 


Thus,  the  layer  strain-stress  relations  (5.47)  may  be  written  as 


k  nk  k 

£  .  =  D  .  .  T  . 

1  1J  J 


(i,  j  =  1,  2,. ..,5) 


(5.73) 


The  layer  stress  components  are  related  to  the  resultants  by 


xk  =  Nkj  (z)nj  (1-1,  2 . 5;  j  =  1 ,  2, ... ,8) 


where  the  matrix  Nv.(z)  is  given  in  Appendix  A. 

*  0 

The  complementary  energy  (5.50)  can  now  be  written  as 


1  L  k  k,+1  k  k 

2  (  S  ^  [  ^ir  <z>Njs  <z>dz3  "r  ", 


(5.74) 


(5  .75) 


B  ( H] »  )  =  ^  ^  rs  nr^ni 9  ^2^s^l5  ^2^  ~  1*  2,. ,.8)  (5,76) 


where  the  result  of  the  thickness  integration  Trs  is  given  in  Appendix  V 
Using  the  above  form  for  the  complementary  energy  density  in  Equa¬ 
tion  (5.44),  one  obtains 


IM 


L  ffe- 


ar 

_  _  _ rs  n 

r  "  3n  "< 


(r,  s  =  1,  2,. ..,8) 


(5.77) 


Taking  the  inverse,  one  obtains  the  resultant-reference  strain 
relations 


nr  -  Eg  (r,  s  -  1,  2,. ..,8) 


(5.78) 


2.  Finite  Element  Formulation 

The  potential  energy  for  the  present  numerical  formulation  is 


TT_  = 


N 

l 

n=l 


[I  [?  Ers  Er  Es  +  ? 


r.  r.  r. 
ij  i  J 


-  Pi  r.]dS 


(r,  s  =  1,  2,... 8;  i,  j  =  1,  2,... 15) 


(5.79) 


where  N  =  number  of  elements,  and  where  the  initial  stress  matrix  N? .  and 

*  J 

the  incremental  load  matrix  Pi  are  defined  in  Appendix  A. 

In  Equation  (5.79),  the  vector  ^  (i  =  1,  15)  refers  to  the  gener¬ 
alized  nodal  displacements  based  on  a  linear  interpolation  scheme  for  the 
three  rectilinear  displacements  and  for  the  two  rotations. 

Thus,  the  reference  surface  strains  can  be  put  in  the  form  (see 
Appendix  A), 

er  *  Brs(V  r,2)rs  (5-80) 


so  that 


N  i-l  * 

=  £  [t  k 
p  n=l  £ 


r  r 
rs  r  s 


PS  rs]n  (r>  5 


1,  2 . 15) 


(5.81) 
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where  the  stiffness  and  load  matrices  k  .  and  p„  in  the  local  frame  are 

rs_  _ s _ _ _ 

k.  =  //  [E  B  .  B.  +  N®  ]dn,  dru 
tq  L  rs  rt  sq  tq  1  2 

and 

Pr  =  II  Pr  dPi  drl2  (5.82) 

Sn 

For  the  rotational ly  symmetric  problem,  the  global  coordinate  frame 
is  that  of  the  lines  of  curvature  coordinates.  Thus,  if  the  transforma¬ 
tion  matrix  from  the  local  to  the  global  frame  is  Trs,  then 


r.  =  T. .  q • 
i  1J 


(5.83) 


where  q.  and  p.  refer  to  the  generalized  nodal  displacements  and  nodal 
T  J 

loads  in  the  global  (lines  of  curvature)  coordinate  system.  Thus, 


?,  C?  krs  <r  "s  -  »r  "A 


(5.84) 


where  the  element  stiffness  and  load  matrices  in  the  global  frame  are 


krs  *  Tri  kij  Tjs 


Pr  ~  "^rs  Ps 


(5.85) 


The  numbering  scheme  for  the  generalized  nodal  displacements  is 
defined  by 
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(5.1 


where  P  is  the  first  node  of  the  triangular  element. 


After  solving  the  variational  equation  671^  =  0  for  the  un¬ 


restrained  generalized  nodal  displacements,  the  following  calcu¬ 
lations  can  easily  be  performed: 

(a)  local  displacements: 


r  =  T  q 
s  sr  ^r 


(5. 


(b)  reference  surface  strains: 
er  =  Brs(rV  r,2) rs 


(5. 


(c)  stress  resultants: 


nr  *  Ers  W  r,2> 


(5. 


(d)  layer  stresses: 

Ti  =  Nij (z)  nj n2)  (5‘ 

(e)  update  the  geometry:  in  the  reference  configuration, 
the  nodal  position  vector  is  given  by  Equation  (3.1), 

r  (P)  =  xi  (P)  e^^  (5. 

thus,  the  new  nodal  position  becomes: 

R  (P)  =  r  (P)  +  [u  (P)  *  e.  1  e. 


(5. 


3 .  Equilibrium  Check 

As  outlined  in  Section  II,  the  equilibrium  check  is  per¬ 
formed  at  each  increment  by  retaining  the  correction  term,  e*, 
defined  by  Equation  (2.20)  in  the  potential  energy  formulation. 
Thus,  Equation  (5.79)  is  modified  as  follows: 


Tf 


P 


N  ft 

SI  I  [h  E  e  e  +  H  N.  .  r.  r. 

JJ  1  rs  r  s  iD  1  3 

n=i  sn 

+  ff  K  £r  -  rr^dS 
s 


p.  rjds 


(5.93) 


where 

n°  =  column (n®x,  n°2>  nj2#  m^,  m°2,  mj2,  r°,  r°) 


and  p°  is  the  load  vector  at  the  previous  increment.  Using 
Equation  (5.80)  for  the  reference  surface  strains,  the  above 
functionsal  becomes 


N  r  A 

"p  =  S  L*  krs  rr 

r  n=l  L 


ys  /\  -i 

r  -  p  r  +  R  r 
s  *s  s  s  sjn 


(5.94) 


where  the  residual  load  vector,  Rg,  is 


Rs  - 


-ff  (nrBrs-^)dS 


(5.95) 


n  ' 


Thus,  the  element  stiffness  equation  becomes 


/v 


k  r  =  p  +  R„ 
rs  s  rr  r 


(5.96) 
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where  is  known  from  the  previous  step  and  p^  is  the  current 
incremental  load. 

Next,  Equation  (5.96)  is  evaluated  for  each  element  and  by 
applying  interelement  compatibility,  these  can  be  assembled  “into--., 
a  system  of  linear  incremental  equilibrium  equations  for  the 
entire  structure, 

kN  %  =  qn  +  *n  <5-97 


whe  re 

KN  =  incremental  stiffness  at  the  Nth  step, 
qN  =  incremental  nodal  displacements, 

Qn  =  current  incremental  load,  and 
=  residual  load. 

If  the  N-l  state  is  exactly  in  equilibrium,  then  the  re¬ 
sidual  load,  R^  is  zero. 

For  the  first  incremental  step  (linear  problem)  R^  =  0 ,  so 

that 

q1  =  (K1)"1  Q1  (5.98) 

After  this  and  every  succeeding  step,  the  geometrical  con¬ 
figuration  is  updated  and  the  total  stresses  and  strains  are  cal¬ 
culated  by  adding  all  incremental  contributions,  and  a  new  stiff¬ 
ness  matrix  is  formed. 
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As  shown  in  Figure  11,  q1  ?  q*,  where  q*  represents  the  true 
solution.  Thus,  the  residual  load,  R2,.  is  calculated  using  the 
solution  q^,  and  at  the  second  step  the  following  system  of 
linear  equations  is  solved: 

K2  q2  =  Q2  +  r2  (5‘" 

The  quantity,  R2,  represents  the  unbalance  in  nodal  equilib¬ 
rium,  introduced  by  the  linearization  process.  This  process  is 

/ 

illustrated  by  Figure  11. 

^r 


Displacement 

Figure  11.  Incremental  Process 
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D. 


Numerical  Results 


First,  a  clamped  square  plate  under  uniform  pressure  is 
considered.  The  variation  of  the  central  deflection  with  respect 
to  the  pressure  load  is  shown  in  Figure  12.  The  exact  solution 
indicated  in  this  figure  does  not  take  shear  deformation  into 
account. 

The  next  example  is  that  of  the  inflation  of  a  six-ply 
aircraft  tire  [7].  Uniform  thickness  is  assumed  with  the  fol¬ 
lowing  material  characteristics: 

Cord  Properties: 

E  =  1.56  x  105  psi, 

vc  =  0.5  psi,  and 

d  (cord  diameter)  =  0.031  in. 

Rubber  Properties: 

Er  =  450  psi,  and 
v  r  =  0.49. 

Tire  Construction  Parameters: 
n  (number  of  plies)  =  6, 

Green  angle  =  57  degrees, 

bead  radius  =  9.15  in., 

number  of  cords  per  inch  =  26,  and 

ply  thickness  =  0.043  in. 


EIS& 


The  undeformed  meridian  profile  is  depicted  by  Figure  13 
and  also  described  in  a  cartesion  reference  frame  in  Table  1. 

The  initial  cord  angle  varies  with  the  meridian  position 
according  to  the  classical  lift  equation  of  bias  tire  construc¬ 
tions,  whose  results  are  shown  in  Figure  14. 

In  the  numerical  solution,  rotational  symmetry  is  employed 
by  specifying  zero  rotation  with  respect  to  the  meridian  and  re¬ 
straining  the  motion  in  the  parallel  direction  of  a  discretized 
strip  consisting  of  80  elements  and  62  nodes,  as  shown  in 
Figure  15. 

The  variation  of  the  crown  displacement  with  respect  to 
inflation  pressure  is  depicted  in  Figure  16,  showing  a  substan¬ 
tial  deviation  from  that  of  the  exact  numerical  solution  of  [7], 
which  does  not  consider  shear  deformation.  Therefore,  further 
verification  is  required  to  evaluate  the  presented  numerical 
results . 
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Figure  12.  Center  Deflection  versus  Pressure  Loading  for  a 
Square  Clamped  Plate 
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Figure  14.  Cord  Angle  Distribution 


94 


asiia."  ■  ijMM in  i - - 


APPENDIX  A 


ATTACHMENT  TO  NUMERICAL  FORMULATION 

Equation  (5.74): 


a  -  2  sk/A;  b  =  2  skz/B;  c  =  -  [  [sk(z!  -  c’)  +  Sk] 


..k 

Nn  » a 

Nu  ■ » 

wk 

N^2  =  a 

4  =  b 

N33  =  a 

N26  =  b 

N47  =  c 

N58  =  c 

(A.  1 ) 

and  all  others  are  zero. 

Equation  (5.76): 

“i  =  4  sk  <ck+i  -  ck>/A' 

a2  =  2  sk  (ck+l  ’  ck^<AB> 

“3  '  f  sk  *ck+l  '  Ck*/B 

“4  =  {sk  »ck+l  -  ck>/5  '  J  ‘k<ck+l  -  ck>  +  ck<Vl  -  ck)] 

+  Sk<ck+1  -  ck>  *  2  sk  Skt<ck*l  '  ck)/3 
'  ck(ck+l  '  Ck>},/B*  (A.2) 
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From  Equation  (5.79),  one  defines 


“elf 

’cn  ci2  ci3_ 

I** 

*Tn" 

e22 

= 

C12  c22  c23 

t22 

1 

CVJ 

_ 1 

-c13  c23  c33- 

-t12- 

and 


k 

k 

r 

L 

y13 

S11  S12 

t13 

_y23_ 

_S21  S22. 

_t23. 

1/  1/ 

where  the  coefficients  matrices  will  be  denoted  by  C  „  and  S  , 

rs  rs 

Then,  it  follows  that 


L 

E 

k=l 


a. 


rr,s+3 


L 

E 

k=l 


a 


k  rk 

2  Vs 


rr+3,s  rs,r+3 
L 

rr+3,s+3  =  kf1  a3  Crs 


(r,s  =  1,  2,  3) 


and 


V+6,s+6 


L 

E 

k=l 


a. 


(r,s  =  1,2) 


(A. 3) 

(A. 4) 

respectively. 


(A.5) 

(A. 6) 


and  all  others  are  zero. 


Equation  (5.80): 

In  a  local  barycentric  coordinate  system,  the  membrane  displacements 
(u] ,  u2),  the  rotations  (ay  a>2),  and  the  normal  displacement  (w)  are 
linearly  interpolated.  Thus,  the  components  of  the  reference  surface  strain 
tensor  are 

e.  =  F..  8.  (i  =  1,2,. ..*1;  j  =  1,2, ...12)  (A.7) 

l  j  j 


where 


e.  =  column  (e^,  e22,  Y-j 2 * 

kH»  k22’  2  k12’  y13*  y23^ 

(A. 8) 

r— 

II 

OJ 

Li. 

II 

Li- 

F32  1  F33  =  1 

F46  =  1  F5,10  =  1 

F6,7  =  1  F6,9  ”  1 

F7,5  =  1  F7,6  =  nl 

F7,7  =  n2  F7,ll  =  1 

F8,8  =  1  F8,9  =  nl 

F8,10  =  n2  F8 ,1 2  =  1 

(A.9) 

The  vector  8i  is  defined  by 

u-j  ~  a-j  +  $i  n-|  +  B2  n2 

t 

u2  =  b-j  +  63  +  n2 

0)i  =  85  +  .Bg  n1  +  s7  n2 
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w2  =  S8  +  89  nl  +  B10  n2 


w  =  c-j  ■*■  3 1  -j  i^i  3 -j 2  h2 


In  the  local  frame,  the  generalized  nodal  displacement  schedule  is 


U^P)  =  r1 

u1  (Q)  =  r2 

u-j  (R)  =  r3 

u2(P)  -  r4 

u2(Q)  =  r5 

u2(R)  =  r6 

^(P)  =  r7 

<^(Q)  =  r8 

w-|  (R)  =  rg 

^2(p)  =  rio 

w2 { Q )  =  r ]  i 

w2(R)  =  r-|2 

w(P)  *  r13 

w(Q)  »  r14 

w(R)  *  r,5 

(A. 1 1 ) 

The  result  of  the 

interpolation  according  to  Equations 

(A. 10)  and 

is 

■  Ru  rj 

(i  =  1,2,. ..12;  j 

=  1 ,2,. ..15) 

(A. 12) 

Thus,  Equation  (A. 7)  becomes 


ei  =  Fi j  Rjs  rs 


so  that  in  Equation  (5 -80 ) » 

Brs^l’  r,2^  =  Frk^nr  VRks 


(A. 13) 


(A. 14) 
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Equation  (5.79): 


The  initial  stress  matrix  in  Equation  (5.79)  is  calculated  from 


c  =  na$  w,a  w,(3 


(A. 15) 


so  that 


Rrp  npk  Rks 


(r,s  =  1,2,. ..15;  p,k  =  1,2,. ..12)  (A. 16) 


wher< 


nl°l,ll 


'11,12 


'12,11 


'12,12 


(A. 17) 


and  all  others  are  zero. 

The  load  matrix  in  Equation  (5.79)  is  derived  as  follows.  For 
constant  pressure  loading,  consider 


pwfa^  n2)  =  P(ci  +  8U  +  612  n2) 


(A. 18) 


In  a  barycentric  coordinate  system,  it  is  sufficient  to  consider 


pw(nr  n2)  = 


(A. 19) 
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where  c-j  is  determined  from  a  linear  interpolation  in  the  form 


C1  =  !11  r13  +  *12  r14  +  *13  r15 


(A. 20) 


Thus,  the  components  of  the  pressure  load  matrix  are 


p13  =  p  *11 


p14  =  p  *12 
p15  =  p  *13 


(A. 21) 


For  the  derivation  of  the  inertial  load  matrix,  begin  with  Equation 
(3.43), 


f  =  [  -  x  /  y(C3)d  Cje 
a  (h)  3  3a 


(A. 22) 


The  mass  density  in  the  above  equation  is  calculated  for  each 
layer  from  the  cord  and  rubber  mass  densities  as 


\ =  y  Yk +  (1  - 


(A. 23) 


where  r^  is  defined  by  Equation  (4.10).  Thus,  the  total  mass  density  m  is 


m  -  /  y(?3)d  C3  -  ^  ^|<(ck+l  ”  ck^ 
(h)  3  3  k=l  K  K  1 


(A. 24) 


wMch  is  concentrated  at  the  centroid  of  each  element,  so  that 


f  =  -  m  Jr  x„  e 


a  a 


(A. 25) 
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APPENDIX  B 

CLOSED  FORM  PLATE  SOLUTION 


For  comparison  purposes,  an  analytical  solution  is  derived  here  for 
the  response  of  a  moderately  thick, simple  supported  plate  under  sinusoidal 

pressure  loading. 


From  Section  H,  the  field  equations  of  an  isotropic  plate  are  as 


follows. 


Equilibrium  equations: 


moS,0  •  ra  ‘  ° 


r  +  p  =  0 

a, a 


(B.l) 
(B.  2) 


Resultant-displacement  and  rotation  gradient  relations, 


mll  =  Dh,l  +  v  w2,2} 
m22  =  D(m252  +  v 

]2  =  D(1  -  v)(wlf2  +  w2>1 )/2 

=  5/h2  D(1  -  v)(Wj  +  c^) 

r2  =  5/h2  D(1  -  v)(wj2  +  w2) 


m 


(B.  3) 


(B.4) 


where 


D  = 


E  h3 

12(1  -  v2)" 


(B.5) 


Mm** 


a&Li&d ;'i  nri  i  -I  Miifai 


104 


The  solution  of  the  above  set  of  equations  may  easily  be  obtained  by 
the  following  change  in  variables: 

$  *  Wj  2  -  Wg  i  (definition) 

©  =  (d-|  i  +  o)2  2  (definition)  (B.6) 

It  can  be  shown  that  the  governing  equations  are 

D  V2  ©  +  p  =  0 

V2  $  =  10/h2  * 

5/h 2  D( 1  -  v)(V2  w  +  ©)  +  p  =  0  (B.7) 

Consider  now  a  rectangular  plate  of  dimensions  a  x  b,  such  that  the 
edges  are  defined  by  x-j  =  ±  a/2  and  X2  =  ±  b/2.  Let  the  pressure  load 
be  defined  by 

p  *  p  cos  ax-j  cos  8  X2 

a  =  fr/a 

8  =  Tt/b  (B.8) 

In  a  simple  supported  situation,  $  is  identically  zero  and  w  and 
©  have  the  forms 

w  =  wQ  cos  a  x-|  cos  8  Xg 

©  =  eQ  cos  a  x-j  cos  8  x2  (B.9) 
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where 


„o  -  P0/D  (1  +  k)  (aJ  l 

So  =  »o/D 


and 


k  =  (a2  +  B2 ) h2/5 (1  -  v) 

The  rotations  and  shear  strains  then  become 
Pn  « 

U1  =  D(a2  +  Sin  01  X1  C°S  3  X2 

P0  a 

w2  =  D(a2  +  02')"r  C0S  a  X1  510  3  X2 


k  a  P0 

*1  =  '  D(a2  +  FT7 

k  e  p0 

y2  D(a2  +  02) 


sin  a  x-j  cos  $  x2 


cos  a  x-j  sin  6  x2 


(B.10) 


(B.ll) 


so  that  the  k  =  0  situation  corresponds  to  the 
The  resultants  become 

mll  =  P0<“*  +  v  32)/(oiZ  +  62)2  cos  a  X1 

m22  =  pQ(v  a2  +  02)/(a2  +  02)2  cos  a  x] 

m]2  =  -  pQ(l  -  v)a8/(a2  +  02)2  sin  a  x] 


Kirchoff-Love  plate  theory. 


cos  6  x2 
cos  0  x2 
sin  0  x2 
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